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Abstract. We consider an unbounded steady-state flow of viscous fluid over a three-dimensional finite body 
or configuration of bodies. For the purpose of solving this flow problem numerically, we discretize the governing 
equations (Navier- Stokes) on a finite-difference grid. The grid obviously cannot stretch from the body up to infinity, 
because the number of the discrete variables in that case would not be finite. Therefore, prior to the discretization 
we truncate the original unbounded flow domain by introducing some artificial computational boundary at a finite 
distance of the body. Typically, the artificial boundary is introduced in a natural way as the external boundary of 
the domain covered by the grid. 

The flow problem formulated only on the finite computational domain rather than on the original infinite domain 
is clearly subdefinite unless some artificial boundary conditions (ABC’s) are specified at the external computational 
boundary. Similarly, the discretized flow problem is subdefinite (i.e., lacks equations with respect to unknowns) unless 
a special closing procedure is implemented at this artificial boundary. The closing procedure in the discrete case is 
called the ABC’s as well. 

In this paper, we present an innovative approach to constructing highly accurate ABC’s for three-dimensional flow 
computations. The approach extends our previous technique developed for the two-dimensional case; it employs the 
finite-difference counterparts to Calderon’s pseudodifferential boundary projections calculated in the framework of the 
difference potentials method (DPM) by Ryaben’kii. The resulting ABC’s appear spatially nonlocal but particularly 
easy to implement along with the existing solvers. 

The new boundary conditions have been successfully combined with the NASA-developed production code TLNS3D 
and used for the analysis of wing-shaped configurations in subsonic (including incompressible limit) and transonic flow 
regimes. As demonstrated by the computational experiments and comparisons with the standard (local) methods, 
the DPM-based ABC’s allow one to greatly reduce the size of the computational domain while still maintaining high 
accuracy of the numerical solution. Moreover, they may provide for a noticeable increase of the convergence rate of 
multigrid iterations. 

Key words. External flows, infinite-domain problems, artificial boundary conditions, far-field linearization, 
boundary projections, generalized potentials, difference potentials method, auxiliary problem, separation of variables. 

AMS subject classifications. 65N99, 76M25 

1. Introduction. 

1.1. Preliminaries. A standard approach to solving infinite-domain boundary- value prob- 
lems on the computer involves truncation as a first step, prior to the discretization of the contin- 
uous problem and the solution of the resulting discrete system. The truncated problem in both 
continuous and discrete formulations is clearly subdefinite unless supplemented by the proper clos- 
ing procedure at the outer computational boundary. The latter boundary is often called artificial 
emphasizing the fact that it originates from the numerical limitations (the discrete system should 
contain no more than a finite number of variables) rather than original formulation. Typically, the 
artificial boundary is introduced as an external boundary of the finite computational domain (i.e., 
the domain covered by the grid, on which the original system is discretized). The corresponding 
closing procedure at the outer boundary is called the artificial boundary conditions (ABC's). 

In the ideal case, the ABC’s would be specified so that the solution on the truncated domain 
coincides with the corresponding fragment of the original infinite-domain solution. However, in 
spite of the fact that different ABC’s methodologies have been studied extensively over the recent 
two decades, the construction of such ideal (i.e., exact) ABC’s that would at the same time be 
computationally inexpensive, easy to implement, and geometrically universal, still remains a fairly 
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remote possibility. The primary reason for that is that the exact ABC’s are typically nonlocal, 
for steady-state problems in space and for time-dependent problems also in time. The exceptions 
are rare and, as a. rule, restricted to the model one-dimensional examples. From the viewpoint 
of computing, nonlocality may imply cumbersomeness and high cost. Moreover, as the standard 
apparatus for deriving the exact ABC’s involves integral transforms (along the boundary) and pseu- 
dodifferential operators, such boundary conditions can be obtained easily only for the boundaries 
of regular shape. 

On the other hand, highly accurate ABC’s are most demanded in many areas of scientific 
computing because as shown by the different authors both theoretically and computationally, the 
overall accuracy and performance of numerical algorithms, as well as interpretation of the numerical 
results, strongly depend on the proper treatment of external boundaries. The possible range of 
applications for different ABC’s techniques is broad. Aerodynamics, in which external problems 
represent a wide class of important applications, especially when it comes to the analysis of three- 
dimensional configurations, constitutes a fraction of this range. Besides the hydro- and aerodynamic 
problems (external flows, duct flows, boundary layers, free surfaces, etc.), the entire range includes 
the flows in porous media, filtration, magneto-hydrodynamic flows, plasma (e.g., solar wind), the 
problems of solid mechanics (in particular, elasticity and aeroelasticity), and the problems of wave 
propagation (electromagnetic, acoustic, seismic), just to name a few. 

As mentioned above, the other usual requirements of ABC’s, besides minimization of the error 
associated with the domain truncation, are low computational cost, geometric universality (i.e., 
applicability to a. variety of irregular boundaries often encountered in real-life settings), and imple- 
mentation without difficulties, in particular, readiness in combining the ABC’s with the existing 
(interior) solvers. The requirements of this group are typically met by many approximate local 
methods that are considered an alternative to the exact ABC’s as the latter are not attainable 
routinely. However, the basic trend in terms of accuracy remains the following: higher accuracy 
for the boundary procedure requires more of the nonlocal nature of exact ABC’s to be somehow 
taken into account. 

In fact, almost any numerical algorithm for setting the ABC’s can be thought of as a com- 
promise between the two foregoing groups of requirements that in a certain sense contradict one 
another. Shifting the balance towards locality and practical efficacy often implies insufficient ac- 
curacy; shifting it to the other side, towards highly accurate nonlocal techniques, may often yield 
cumbersome and all but impractical algorithms. It is not surprising, therefore, that the treatment 
of external boundaries in modern production computations typically follows the first, local, path. 
In computational fluid dynamics (CFD), for example, only a few ABC’s methodologies out of the 
wide variety proposed to date can be regarded as the commonly used tools. All of them are ei- 
ther based on the essential model simplifications, e.g., local quasi-one-dimensional treatment in the 
vicinity of the artificial boundary, or obtained as a localization of some nonlocal ABC’s. To meet 
the overall accuracy requirements when using such simple boundary procedures, one often has to 
choose the excessively large computational domains. 

A survey of methods for setting the ABC’s in different areas of scientific computing can be 
found in our work [1], as well as in the comprehensive reviews by Givoli [2, 3]. These surveys give 
a comparative assessment of advantages and disadvantages of various global and local techniques, 
and also point out the relations between the global and local methods. 

1.2. Methods and Objectives. This paper continues our work on constructing the ABC’s 
that would combine the advantages relevant to both local and nonlocal approaches. The specific 
area of applications that we are looking at is steady-state external viscous flows. 

Previously, we have developed and implemented in practice the highly accurate ABC’s for two- 
dimensional case (plane geometry). Our approach is based on usage of the Calderon generalized 
potentials and pseudodifferential boundary projections [4] (see also work by Seeley [5]). The po- 
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tentials and projections are actually employed in the modified form proposed by Ryaben’kii; the 
corresponding numerical technique used for calculating the potentials and projections is known as 
the difference potentials method ( DPM ), see work by Ryaben’kii [6, 7, 8] and also the description 
of the method in the book by Mikhlin, et al. [9]. The resulting DPM-based boundary conditions 
appear global. As, however, will be seen, one of the principal advantages that we gain using the 
DPM is that the method allows us to simultaneously meet the high accuracy standards and the 
requirements of geometric universality and easiness in implementation. 

The two-dimensional DPM-based ABC’s have been used along with the multigrid Navier-Stokes 
code FLOMG by Swanson and Turkel [10, 11, 12]. In spite of their nonlocal nature, the new boundary 
conditions readily apply to the boundaries of irregular shape and appear very easy to incorporate 
into the existing solver. In our computations, the DPM-based ABC’s have explicitly outperformed 
the standard local methods from the standpoints of accuracy, convergence rate, and robustness . 
The investigated regimes range from the very low (incompressible limit) to transonic Mach numbers 
and encompass both laminar and turbulent flows. 

The aforementioned two-dimensional constructions and the corresponding numerical results 
have been reported in a series of papers. In work [13], we describe the foundations of the DPM- 
based approach to setting the ABC’s for computation of two-dimensional external viscous flows 
(Navier-Stokes equations). In work [14], we implement this approach along with the code FLOMG 
and present some numerical results for subsonic and transonic laminar flows over single- element 
airfoils. In work [15], we show the results of subsequent numerical experiments and propose an 
approximate treatment of turbulence in the far field. Our work [16] delineates the algorithm for 
solving one- dimensional systems of ordinary difference equations that arise when calculating the 
generalized difference potentials. In work [17], we extend the area of applications for the DPM-based 
ABC’s by analyzing two-dimensional flows that oscillate in time; we also provide some solvability 
results for the linearized thin-layer equations used for constructing the ABC’s. In work [18], we 
present a general survey of the DPM-based methodology as applied to solving external problems 
in CFD, including parallel implementation of the algorithm, combined implementation of nonlocal 
ABC’s with multigrid, and entry-wise interpolation of the matrices of boundary operators with 
respect to the Mach number and the angle of attack. Additionally, in [18] one can find some new 
theoretical results on the computation of generalized potentials, the construction of ABC’s based 
on the direct implementation of boundary projections (thin-layer equations), and some numerical 
results for various airfoil flows: laminar and turbulent, transonic and subsonic, including very low 
Mach numbers. 

The next natural objective after constructing the two-dimensional algorithm is the analysis 
of three-dimensional steady-state flows. This case is undoubtedly the one most demanded by the 
current practice in CFD. In work [19, 20], we outline the basic elements of the DPM-based ABC’s 
for steady-state viscous flows around the wing-shaped configurations and show some preliminary 
numerical results for the subsonic regime. The numerical results of work [20] are obtained with the 
NASA-developed production code TLNS3D by Vatsa, et al. [21]. In work [22] we further develop 
the three-dimensional DPM-based algorithm and present the computational results for transonic 
flows. In all cases (see [20, 22]), the DPM-based ABC’s allow one to greatly reduce the size of the 
computational domain (compared to the standard local boundary conditions) while still maintaining 
high accuracy of the numerical solution. This actually means the overall increase of accuracy due 
to the improved treatment of the artificial boundary; it also implies the substantial economy of the 
computer resources. Moreover, the DPM-based ABC’s may provide for a noticeable speedup (up 
to a factor of three) of the convergence of multigrid iterations. 

Below, we for the first time systematically describe the three-dimensional DPM-based ABC’s 
for calculating viscous flows around the wings. We address the theoretical foundations of the 
approach, present the numerical algorithm with a fair amount of details, and demonstrate the 
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computational results for the different flow regimes, including the low speed flow and the flow with 
the developed separation. The numerical results for the DPM-based ABC’s are compared with 
those obtained with the standard local method. 

The material in the paper is prepared as follows. In Section 2, we formulate the problem, 
describe the specific geometric setting for three space dimensions, provide the foundations for the 
DPM-based ABC’s on the continuous level, and then implement the new algorithm in the finite- 
difference framework. In Section 3, we first briefly summarize the results of our previous numerical 
experiments in two space dimensions and then report on the recent three-dimensional computations 
for various flow regimes. Section 4 contains the conclusions. 

2. External Flow. 

2.1. Formulation of the Problem. We consider an unbounded steady-state flow of viscous 
fluid past a three-dimensional wing. The flow is uniform at infinity. We consider both incompress- 
ible and compressible formulations, in the latter case we assume that the fluid (gas) is thermody- 
namically perfect and that the free stream is subsonic. Moreover, as the fluid is viscous and the 
size of the immersed body (wing) is finite, the flow limit at infinity is a free stream. 

Generally, the near-field flow is governed by the full Navier-Stokes equations. However, in 
many cases (including those studied in this paper, see Section 2.2) the full system can be simplified 
and reduced to the so-called thin-layer equations [23], which do not contain stream wise viscous 
derivatives. In particular, this simplification is done in the code TLNS3D that we are using for our 
numerical tests (Section 3). Moreover, for the most interesting case of turbulent flows the near- 
field numerical algorithm should also involve some turbulence model, we comment on this issue in 
Section 3, which is devoted to numerics. 

2.1.1. Linearization. In the far field (i.e., far enough from the finite immersed body), the 
perturbations of the flow induced by the immersed body are small and we therefore linearize 
the governing thin-layer equations against the constant free-stream background. Introducing the 
Cartesian coordinates ( x , y, z) and assuming (without loss of generality) that the free stream is 
aligned with the positive x direction, we can write the dimensionless linearized equations as 


(2.1a) 
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du 
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Here, u, v , w , and p are the perturbations of the Cartesian velocity components and pressure, 
respectively, and p is the perturbation of density for the compressible case. 

The derivation of equations (2.1) involves two consecutive steps. First, we introduce the di- 
mensionless quantities for the original thin-layer system. To do that, in both incompressible and 
compressible case we scale the velocity projections by the dimensional free-stream value of the 
^-component u$. (As mentioned above, the free-stream y and z velocity components are zero: 
Uo=0, u?o=0.) Moreover, for the incompressible flow we scale the original pressure by no 2 and for 
the compressible flow the pressure is scaled by p 0 uo 2 , the internal energy e is scaled by uo 2 and 
the viscosity p is scaled by p 0 . (Everywhere above, the subscript “0” denotes the corresponding 
“full,” i.e., “nonlinear,” dimensional value.) Finally, we scale the coordinates x, y, and £ by the 
characteristic length X, for example, it may be the root chord or the semi-span of the wing (see 
Section 2.2). For the compressible flow, we also have to use the equation of state of the perfect gas 
to eliminate the internal energy from the original system. 

After the nondimensionalization, we represent each quantity (velocities and pressure for the 
incompressible case and velocities, pressure, and density for the compressible case) as a sum of the 
constant background value (free stream) and the small perturbation and retain in the equations 
only the first-order terms with respect to the perturbations. (Note, in the incompressible case only 
the gradient of the pressure is involved in the original system, therefore the actual value of the 
background constant for the pressure does not matter.) In so doing, we arrive at equations (2.1a), 
(2.1b) for the incompressible flow and equations (2.1a), (2.1c) for the compressible flow. Both in 
(2.1b) and (2.1c) Re is the Reynolds number (in the turbulent case, it is an effective far-field value, 
see [15]); in addition, for the compressible flow (see (2.1c)) Mo = u 0 / ( IPo/po ) ly ^ 2 is the free-stream 
Mach number (always M 0 < 1), Pr is the Prandtl number, and 7 is the ratio of specific heats. 

Note, for the incompressible case it is clear that the differential equations for the small per- 
turbations are linear. For the compressible case, however, this fact may require some additional 
justification, see Section 2.1.2. 

The system (2.1) describes the flow in the far field. In both incompressible and compressible 
cases, it is supplemented by the boundary condition 


( 2 . 2 ) 


u — *■ o, as r = (x 2 + y 2 + z 2 ) 1 / 2 — * +oo, 
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which means that all the perturbations vanish at infinity, or equivalently, the flow approaches the 
free stream. 

Let us mention, that the matrices C, D , E, F, H of ( 2 . 1 b) are symmetric, whereas the matrices 
C , D, E, F , H of (2.1c) are not. As the symmetric form of the matrices may sometimes be more 
convenient for the analysis and also more suitable for the numerical calculations (especially when 
the Mach number Mo is low), one can use the transformation proposed by Abarbanel and Gottlieb 
in work [ 24 ] to simultaneously reduce C,D,E,F,H of ( 2 . 1 c) to the symmetric (and some of the 
matrices to the diagonal) form. Specifically, introducing the non-degenerate matrix S, 
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we have instead of (2.1c) 
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The symmetric form C, D , j E, JF, IT, J (see (2.4)) of the system matrices appears useful when 
analyzing the incompressible limit Mo — > +0. A detailed study of the compressible Euler and 
Navier-Stokes equations as the Mach number approaches zero can be found in work [25], as well as 
in work [26, 27]. 

2.1.2. Asymptotic Methods — Linear vs. Nonlinear. In connection with the linearized 
model proposed above for the far field, especially as it regards the compressible case, we will discuss 
here one group of the ABC’s methods that are widely used in computations. These methods employ 
some asymptotic form of the far-field solution for closing the system of equations to be solved inside 
the computational domain. Typically, this approximate asymptotic form can be obtained as a few 
leading terms of the series (or asymptotic series) that represents the solution in the far field. The 
corresponding ABC’s most often appear local. 

The idea of the this type was employed, for example, by Sa and Chang in work [28] to set 
the ABC’s for vorticity when integrating the incompressible Navier-Stokes equations around a 
cylinder. Burkhart in [29] and Burkhart et al. in [30] derive an asymptotic expansion for the finite- 
difference fundamental solution of the three-dimensional Laplace operator on a Cartesian grid and 
then use a few leading terms of this expansion to set the ABC’s for an external flow problem that 
is solved within the full-potential framework. Wubs, Boerstoel, and Van der Wees in [31] use a 
Fourier representation of the far-field solution to the two-dimensional Laplace equation to calculate 
a potential flow around an airfoil. The ABC’s [31] are again derived from the first few leading 
terms of the expansion; as the artificial boundary approaches the airfoil, more terms are required 
to maintain the accuracy. Finally, the so-called point-vortex model, which has been proposed by 
Thomas and Salas in an earlier work [32] and which is extensively used in the today’s CFD is also 
based on the idea of asymptotics. Specifically, the first leading term of the far-field expansion for 
the linearized flow potential is used to calculate the velocity projections at the external boundary 
when computing the two-dimensional compressible flows. This leading term is proportional to the 
circulation of the flow. 

The asymptotic methods may often require the explicit knowledge of the coefficients that multi- 
ply the corresponding terms in the expansion (the ones that are used in the far-field representation). 
In CFD, these coefficients are typically obtained through the boundary conditions on the surface 
of the immersed body. For example, the value of circulation for the point-vortex model [32] is 
proportional to the lift, which is calculated by integrating the pressure along the surface. There is 
also another way of using the asymptotics for setting the far-field ABC’s. It has been proposed by 
Ba.yliss and Turkel in work [33, 34, 35] and by Bayliss, Gunzburger, and Turkel in work [36] and 
does not require the explicit knowledge of the coefficients. Instead, the authors of [33, 34, 35, 36] 
develop a set of special local differential relations that identically cancel out the prescribed number 
of leading terms in the corresponding series; these relations can obviously serve as the ABC’s. 
However, they are typically of a high order even if the order of the original equation (system) is 
low. 

As a rule, the asymptotic ABC’s methods are derived on the basis of the linear (or linearized) 
equations (apparently because it is easier to study the convergence of the corresponding series). In 
certain cases, however, one takes into account the nonlinear corrections as well. For example, when 
analyzing the transonic limit Mq — > 1 for the small perturbations of the velocity potential in the 
flow of compressible gas, some second-order terms should formally be retained in the differential 
equation along with the first-order terms (see, e.g., the book by Cole and Cook [37]). This leads 
to the nonlinear Karman-Guderley equation rather than to the linear Prandtl-Glauert equation 
(the latter is valid for smaller Mo). For two-dimensional external flows (e.g., flows around airfoils) 
described by the Karman-Guderley equation, it turns out that the nonlinear corrections to the 
leading linear lift-based term ~ —T0/2ir in the far-field expansion of the potential (T is the flow 
circulation, 8 is the polar angle) contain the terms proportional to log r/r (r is the polar radius), 
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which formally decay slower than the next linear term ~ 1/r as r — > +oo. This circumstance, 
in particular, gave reasons to Drela in work [38] and Giles and Drela in work [39] to include the 
nonlinear correction terms in their simplified far-held potential model for the compressible airfoil 
calculations. (Note, the entire series that represents the behavior at infinity of the potential function 
of a two-dimensional subsonic compressible flow has been accurately constructed by Ludford in work 
[40] using the hodograph plane techniques.) 

Our two-dimensional DPM-based approach of [13], however, uses the full flow system; we never 
introduce the potential and we always consider only the linearized far-held how. The accuracy and 
performance of the corresponding nonlocal ABC’s are demonstrated by the numerical experiments, 
see our work [14, 15, 18]. These accuracy and performance are typically better than those of the 
standard methods. However, we should say that the investigated Mach numbers have never come 
really close to the transonic limit, we have always run our calculations in the range Mo ^0.8. 
Generally, retaining the second-order nonlinear terms in the compressible far-held model for two- 
dimensional hows is most relevant to the case of Mach numbers close to one, Mo — » 1, whereas 

the linear theory works best for 6 <C (l — M 0 2 ) (see [37]), here 6 can be regarded as, e.g., the 
airfoil thickness. 

The situation with the compressible far-held expansion for three space dimensions is entirely 
different. Let us consider here the Karman-Guderley equation (see [37] ) 

d 2 <£ d 2 4> 7 + 1 d<l> 0 2 <?s> 

dx 2 dy 2 dz 2 K dx dx 2 ' 

In equation (2.5), <j) is the perturbation of the full potential $ of the flow around a thin three- 
dimensional wing so that 


(2.6a) 


J_0* = 1 ,2/3^ i_^> 

uq dx dx ’ uq dy dy ’ 


Uq dz dz ’ 


y = $ 1/3 y, y = 6 1/3 y, 


6 is the wing thickness (6 — - +0 along with Mq — * 1 in the transonic limit), 


(2.6b) 


. 1 - Mq 2 

' £ 2/3 


is the parameter of transonic similarity (the true linear theory corresponds to big values of K , see, 
e.g., [41]), the additional coordinate transformation is given by 


(2.6c) y — VKy, z = \[Kz , 

7 is still the ratio of specific heats, and the free stream is again aligned with the positive x direction. 

The far-held expansion for 0 in the linear theory, i.e., when the right-hand side (RHS) of 
equation (2.5) is omitted, starts with the horseshoe vortex (see, e.g., [42]) 


(2.7a) 


_ y ( _ (1 + cos 0) cosy? 

01 ~ f + p\ i+ r) “ f sin 0 


where the spherical coordinates are introduced as 
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x = r cos 6 , y = r sin 6 cos <p, z — f sin 6 sin <p. 


Expression (2.7a) obviously has a singularity in the wake, i.e., along the line 0 = 0. Clearly, the 
source term of order 1/r is not present in the far-held expansion because the surface of the wing is 
assumed closed. Therefore, the next term in the linear expansion should be proportional to 1/f 2 . 
We consider its general form 4> 2 = X/,™ y?), l = 1, m = -1, 0, 1, where the spherical 


functions Y \ m (0,(p) are given by Y t m (6, <p) = V™ (cos 9)exp(im<p), and = (1 - p 2 ) m l 2 v\ m \p), 

1 d) 

to < l, here Vi(p) = ~ 1)* are the Legendre polynomials. Using the real representation, 

we obtain the general form of the second term as 


(2.7b) 


, x , y z 

(f>2 ~ a TZ + b-^r + = 


cos 6 . sin 0 cos p sin 6 sin p 

+ b + c — 


where a, 6, and c are some arbitrary constants. 

To obtain the nonlinear corrections due to the RHS of equation (2.5), one can substitute the 
linear terms (2.7) into this RHS and solve the resulting Poisson equation. For the purpose of 
simple demonstration we will do that separately for <fi i of (2.7a) and (j ) 2 of (2.7b) although the 
similar procedure can be carried out for any weighted sum of fa and fa. 

Substitution of (2.7a) into the RHS of equation (2.5) yields 


/0 1 d f 1 d /. n d(j)\ 1 d 2 cf) 7 + 1 3 cos 0 sin 2 0cos 2 p 

(2-8) ^2^ ( ^ TTT ) T — T“ T 7T Sin0 — + To . 2/} T~T = Y T “• 

r z or \ or J r 2 sin 9 86 \ d9 ) r 2 snr 6 dp 2 K r 5 

Note, the singularity of the potential (2.7a) in the wake (along 0 = 0) vanishes with differentiation. 
The RHS of equation (2.8) can now be expanded with respect to the spherical functions; the 
corresponding finite Fourier series is, in fact, given by 


(2.9) cos 6 sin 2 6 cos 2 <p = - ^Y 3 O (0>¥>) + ^Y£(0,<p). 

Taking into account that the spherical functions Yfa are actually the eigenfunctions of the Beltrami 
operator on the sphere (two last terms on the left-hand side of equation (2.8)), we can separate 
the variables in equation (2.8) and reduce it to the finite family of one- dimensional equations with 
respect to the Fourier transformation fa, m = fa 


( 2 . 10 ) 


cl 2 4> 2 d(j) 1(1 + 1)- Ai 

dr 2 f dr r 2 f k + 2 5 


/= 1,3, k = 3. 


The constants A\ in (2.10) are, of course, inverse proportional to the transonic similarity parameter 
Ii of (2.6b), they also involve the coefficients of the expansion (2.9). The homogeneous counterpart 
to equation (2.10) has two linearly independent solutions, <f>i(f) ~ f l and fai(f) — f~ l ~ 1 . The 
solution far) to the nonhomogeneous equation (2.10) can therefore be found in the form 


(2.11a) </>(f) = cj(f)4>i(f) + cn(f)4>u(f) = ci(f)f l + c n (f)r 1 \ 

where the functions cj(r) and cu(f) satisfy the system 
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(2.11b) 


4>i(r) 

d<f>i(r) 


fal(r) 


df dr 

Solving (2.11b) with respect to c/(f) and yields 


d 

' C/(f) ■ 


0 

df 

•-H 

. 1 


Ai 

^/c+2 


( 2 - 12 ) c *(0 ~ c "( f ) ~ 

provided that / + k / 0 and — l + & zfz 1, respectively. As both of the latter conditions are met for / 
and k from (2.10), we substitute the expressions (2.12) into equality (2.11a) and finally obtain the 
nonlinear correction due to the horseshoe potential of (2.7a) as 


(2.13a) (f> lNL ~ 

The same type of derivation can be performed for the doublet potential of (2.7b). Substituting 
<f ) 2 of (2.7b) into the RHS of equation (2.5) one obtains the expression proportional to f“ 7 instead 
of f" 5 in the RHS of equation (2.8). The expansion analogous to (2.9) will now contain for 
/ = 1, 3, 5, and after the separation of variables the equation (2.10) will also change accordingly, 
we will have k = 5 instead of k = 3 and add l = 5 to the set of wavenumbers. As a result, the 
nonlinear correction due to the potential 4 > 2 of (2.7b) can be shown to have the form 

(2.13b) 4> 2nl ~ 1 

We see that the nonlinear correction i NL of (2.13a) decays at infinity two orders of magnitude 
faster than the term <j>\ of (2.7a.) that it originates from. Analogously, the nonlinear correction 
<j> 2 nl of (2.13b) decays at infinity three orders of magnitude faster than the corresponding term 
cf >2 of (2.7b). We therefore conclude, that unlike the two-dimensional case, the transonic nonlinear 
corrections are not required when analyzing the far field for three space dimensions. This conclusion 
basically coincides with the results of [37] saying that the far field around a thin three-dimensional 
finite-span wing is essentioally linear. In other words, we have shown that the small perturbations of 
the velocity potential of a three-dimensional compressible flow are described by the linear formulas 
in the far field even when Mo — > 1. This justifies the far-field linearization of Section 2.1.1. 

2.1.3. Outline of the Algorithm. Having introduced the linearized model (2.1), (2.2) for 
the far filed (see Section 2.1.1) we now obtain the combined problem: nonlinear inside the finite 
computational domain and linear on its infinite exterior. The nonlinear and linear parts of the 
problem are, of course, not independent. The interior and exterior solutions should match at the 
artificial boundary and consequently, the combined problem must be solved as a whole. Therefore, 
at the first glance the new problem is no easier than the original one from the standpoint of 
solving it numerically because it is still formulated on the unbounded domain. However, using the 
methodology of Calderon’s projections and the DPM [6, 7, 8], the exterior linear problem can be 
effectively reduced to a certain nonlocal relation formulated on the artificial boundary. The latter 
relation can serve as the desired ABC’s. 

More precisely, we introduce a special space of (vector-)functions at the artificial boundary, 
it will be called the space of clear traces [6, 7, 8]. We also define a (pseudodifferential) projec- 
tion operator that maps the space of clear traces onto itself. This operator will be analogous to 
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Calderon’s boundary projections. Under certain conditions, one can show that the element of the 
space of clear traces is actually a trace of some solution to the problem (2.1), (2.2) on the exterior 
domain if and only if this element belongs to the image of the aforementioned projection operator. 
The latter condition can, in fact, be written in the form of Calderon’s boundary equation or the 
so-called boundary equation with projection (BEP) [6, 7, 8]. Its solutions will provide us with the 
complete boundary classification (in terms of the appropriate traces) of all those and only those 
tPs (see (2.1b) and (2.1c)) that solve (2.1), (2.2) outside the computational domain. 

As we intend to set the ABC’s for the discretized flow problem, the foregoing boundary classifi- 
cation of the exterior solutions will also be obtained in the discrete framework using the concept of 
finite-difference clear traces and finite-difference counterparts to Calderon’s projections and gener- 
alized potentials [6, 7, 8]. Once we are able to calculate the image of Calderon’s projection (i.e., the 
result of application of this operator to every given input), we can actually set the ABC’s in a few 
different ways. Earlier (see [13, 14, 15]) we have been solving the corresponding BEP using some 
variational approach. Below, we follow the different path, namely, we implement the boundary 
projections directly as proposed first in [18]. In fact, applying the Calderon operator we update 
the missing boundary values on every cycle of the iteration procedure that is employed inside the 
computational domain. 

2.2. Geometric Issues and the Basics of the Discrete Algorithm. The specific config- 
uration of domains that we will be dealing with hereafter is shown in Figure 2.1. 



Fig. 2.1. Schematic geometric setup; the wing on the left is enlarged. 

The actual structure displayed in this figure is the well-known test wing ONERA M6 (the 
wingtip is blunted, it is in the “hidden” area on the figure). The wing stretches span- wise along 
the Cartesian axis £ and is assumed symmetric with respect to the plane z = 0. The uniform at 
infinity fluid flow is coming along the positive x direction, which together with the symmetry of 
the wing implies the symmetry of the entire flow pattern with respect to z = 0. 
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The latter symmetry, in turn, means that the components u‘ of u, i = 1, . . .,4 for (2.1b) and 
i = 1, . . .,5 for (2.1c) (note, the x projection of the velocity vector u = u l for i = 2) satisfy the 
following set of equalities for i ^ 4: 


(2.14a) 


u l (-,\z\) = 


du\ du\ du\ du\ . .. du\ . .. du\ . .. 

dx ^ dy { ' ’ 


C'» l* e 'l) ^„.2 ("’ l z D’ ^ r 2 Y I- 2 !) — Y N)> Y) ” mm~Y 


dy 2KUU dy 2 

which, in particular, implies 


d;? 2 


dydz 


dydz 


{2.14b) ^-(q0) = 0, J^(-,0) = 0, fo r* #4, 

and the following set of equalities for i — 4: 


(2.15a) 


«*YM) = -«*’(•, -kl) 


du‘ . . 


du x * 0tt*‘ . . 

' dx^’ ^ dy^ } 


d 2 ^ 


d 2 u l 


d 2 u l 


du { 

d 2 u l 


dy 2 

which yields 




dz 2 


dz 2 




d 

d 2 u l , | | 




(2.15b) 


«'M) = o, ir ( ' 0) = 0 ' §Y 0) = 0, 

9V (',0) = 0, ^(,0) = 0, fori = 4. 


dy 2 


Relations (2.14), (2.15) will be used in Section 2.4 when constructing the discretization for the 
linearized system. 

The flow equations are integrated numerically on a curvilinear grid generated around the wing. 
The grid shown in Figure 2.1 is a one-block C-0 type grid; it this paper we will use the grids of this 
type only. The surface designated T on Figure 2.1 is actually the external set of nodes of the C-0 
grid, i.e., the artificial boundary. Henceforth, we will also need the notation D{ n for the interior of 
T, i.e., for the finite computational domain, and the notation D ex for the infinite exterior of T. 

The curve Ti C D ex (see Figure 2.1) actually represents the set of ghost nodes (or centers of the 
ghost cells for the case of finite- volume discretization), it can also be thought of as the outermost 
set of nodes of the original C-0 grid; the surface T then becomes the penultimate set of nodes. We 
will further assume that the linearization (2.1) is valid in D ex , i.e., outside F, so that T x belongs 
already to the linear zone. The actual admissible size of Di n such that the perturbations can be 
considered sufficiently small and therefore the assumption of the linearity in D ex would hold, is, of 
course, unknown ahead of time. We verify the validity of linearization in D ex a posteriori, through 
the series of numerical experiments, see Section 3. 
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Clearly, when the stencil of the scheme used inside Di n is applied to any node from T, it 
generally requires some ghost cell data. Note, for the second-order central difference schemes 
(like those employed in the code FLOMG, see [10, 11, 12] and TLNS3D, see [21]) the consideration 
of only one layer of ghost cells is sufficient, but the case when the stencil is more extensive 
and consequently, more ghost cells are required, can be treated similarly. Unless the missing 
ghost cell data are provided, i.e., obtained with the help of the ABC’s, the discrete system solved 
inside the computational domain appears subdefinite, in other words, it has less equations than 
it has unknowns. As mentioned in Section 2.1.3, the apparatus of the DPM [6, 7, 8] gives us the 
complete boundary characterization of the traces of exterior linear solutions. Since the linear zone 
D ex stretches from F to infinity and contains Ifi, the following approach appears most natural 
for setting the ABC’s. First, the data provided from inside D{ n are subjected to the projection 
operation. The resulting projection will by definition admit a complement on D ex that solves 
(2.1), (2.2). The latter complement can be calculated in the form of a generalized potential and 
considered on Ti. Altogether, this procedure yields the missing relations between the values of the 
solution on T and T\. In other words, it provides for a desired closure to the discrete system solved 
inside D { n , or the ABC’s. Typically, the solution algorithm inside D{ n involves some pseudo-time 
iterations (see Section 3); then, the foregoing closing procedure is applied on every iteration cycle, 
more precisely, every time the ghost cells need to be updated in order to advance the next time 
step. 

We now proceed to the description of the generalized potentials and boundary projections, as 
well as their finite-difference counterparts, that are required for setting the DPM-based ABC’s. 
Note, if the potentials and projections are calculated for exactly the operator L of (2.1a) that 
operates on the functions u defined on the entire inffinite domain D ex and satisfying boundary 
conditions (2.2), then the corresponding BEP appears equivalent to the original linear problem 
(2.1), (2.2) (see [6, 7, 8]). We, however, will have to introduce some simplifications and to carry 
out the DPM-based procedure for a certain approximation of problem (2.1), (2.2) (see Sections 2.3 
and 2.4) rather than for this problem itself. Nonetheless, the corresponding approximate solution 
can be made as close to the original one as initially prescribed (see [18] and below for more detail). 
Therefore, within the accuracy of far-field linearization the resulting ABC’s can be made as close 
to the exact ones as desired. 

2.3. Foundations of the DPM-based ABC’s. Here, we will first formulate and solve the 
so-called auxiliary problem (AP) for the inhomogeneous version of system (2.1) with boundary 
conditions (2.2). This AP will be the central element of our construction of Calderon’s generalized 
potentials and boundary projections. In fact, the solution of the AP can be thought of as a 
substitute for the convolution with the fundamental solution in classical potential theory. 

2*3.1. Infinite-Domain AP. Let us consider a compactly supported vector-function {/*} = 
f = f(x,y 7 z), supp / C Di n 5 of dimension four (i = 1,...,4 for (2.1b)) or five (i = 1,...,5 for 
(2.1c)); in the meantime, we do not specify the concrete form of f. This function f will be the 
right-hand side for the AP. The AP is initially formulated on the entire space R 3 ; namely, we will 
be looking for a solution u of system 


(2.16) Lu = / 

that meets boundary condition (2.2) at infinity. 

Note, when discussing regular solutions below we assume, if necessary, that the functions 
involved can be represented in the form of their Fourier integrals. 

Proposition 2.1. Let f be a compactly supported distribution , f £ V'(R 3 ), supp f C D{ n . 
Then , system (2.16), (2.1b) is solvable in the Schwartz space V f (R 3 ). 



14 


S. V. TSYNKOV 


To justify proposition 2.1, we use the standard methodology based on application of the Fourier 
transform over the entire R 3 (see, e.g., [43]). Clearly, to obtain the solvability of (2.16) in V'(R 3 ) 
it is sufficient to make sure that the inverse symbol of the operator L (see (2.1)) belongs locally to 
Li(R 3 ). Denoting the dual (Fourier) variables to ( x , y, z) by (£, rj, £), we can write the symbol Q 
of L as 


(2.17) Q = i£C + ir/D + i(E - tj 2 F - ( 2 H - r](J. 

Then, the entrees qjk = qkj, j = 1, . . .,4, k = 1, . . .,4, of the inverse symbol Q _1 for the incom- 
pressible case (2.1b) are given by 


(2.18) 


qn = 


+ 


Re 


q 22 


v 2 + C 2 




i n 2 +< 2 \ 

+ Re J 




qi2 = — o~, 
g 2 

-(n 

II 

CO 

<723 = 7 

(*« 

+ 1,2 1 <2 ' 

' Re 

e + c 2 

)t> 2 ’ 


1 v 2 +( 2 ' 

"T" Re 

w 


—it) 


q-24 

<734 


q 44 



-*c 

q \4 = 


-K 



j? 

-v( 


(•■f + 

w 

e+v 2 


(i(+^ 

w 


where g 2 = T; 2 + rf + ( 2 . From (2.18) one can see that Q -1 has only one real singularity, it is 
located in the origin, (£, r/, () — (0, 0, 0). Clearly, all difc, k = 1 ,. . .,4, are absolutely integrable 
near the origin and consequently, pu- € L l ° c (R 3 ) for k = 1, . . .,4. As for the other qjk, we introc a 
the polar coordinates: £ = geos 6, r/ = psin0cos<p, ( = psin#sin<^, 0 < 6 < ir, 0 < <p < 2ir. d 
notice that for sufficiently small p’s 


(2.19) 


q2 < 2 a , Q 2 

< cos 0 + 


sin 4 0, 0 < 6 < n. 


Re 2 - Re 2 

For j = 2, . . . , 4 and k = j,..., 4, estimate (2.19) immediately yields 


\qjk\ < const 


const 




Q 2 q ^cos 2 6 + -^2 sin 4 6 s ) 


1/2 


Re 

< const— 


and we therefore conclude that qju £ L l { c (R 3 ) for all j = 1,...,4, k = 1,...,4. Thus, we have 
shown that proposition 2.1 does hold for the incompressible case (see (2.1b)). Note, similar proof 
for the two-dimensional compressible case can be found in [17]. 

Proposition 2.2. Let f £ L\{R 3 ). Then , equation (2.16), (2.1b) may have no more than one 
regular solution u £ V\R 3 ) that vanishes at infinity, i.e. satisfies (2.2). 

Indeed, the Fourier transformation f of the RHS is continuous on R 3 in this case. The regular 
solution u to (2.16), (2.2) is given by the inverse Fourier transform: u — (Q~ l • Since Q~ r f 
has only one real singularity (in the origin), then any other solution can differ from u by no more 
than an inverse Fourier transformation of a distribution concentrated in the origin. The latter 
can be nothing but a sum of ^-functions and their derivatives, which correspond to polynomials 
after the Fourier transform. Therefore, proposition 2.2, which is actually a statement of conditional 
uniqueness, has been justified. Note that we have been able to establish uniqueness so easily because 
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the inverse symbol Q 1 has only one isolated real singular point. This, for example, would not be 
the case for the Euler equations, which can be obtained by formally letting Re~ x = 0. 

Proposition 2.3. Let f be compactly supported and f G L 2 (R 3 ). Then , equation (2.16), 
(2.1b) has a solution u = u 1 + u 11 , where u 1 is infinitely smooth on R 3 , i.e., u 1 G C°°(R 3 ), and 
satisfies boundary condition (2.2), and u 11 G L 2 {R 3 ). Moreover , for any e > 0, one cam always 
choose the representation u = u 1 + u 11 so that lltt 


II 


L 2 (R z 


< €. 


Consider a partition of unity l = go + go, where both (scalar) functions go and <75 are infinitely 
smooth on R 3 , g 0 = 1 on a ball Ur 0 centered in the origin with the fixed radius R 0 , and g 0 ~ 0 
outside a bigger ball f/# 0+Ai with the radius R 0 T //, > 0. Clearly, as f G L 2 (R 3 ) and f is 

compactly supported, then / G Li(R 3 ). The solution u is given by the inverse Fourier transform: 
U= (Q-'fY = (Q _1 5o/)~ + (Q _1 </o/) • The first term of the foregoing sum, u l = ", 

is obviously an inverse Fourier transformation of a function from Li(R 3 ). Moreover, for any poly- 
nomial P = P(£, g, C) : PQ-'gof € ii(-R 3 ). Therefore, u 1 G C°°(R 3 ) and u 1 meets boundary 
condition (2.2). The second term, u 11 = _i ffo iy , is an inverse Fourier transformation of a 

function from L 2 (R 3 ) because f G L 2 (R 3 ) and Q _1 is bounded when g — > +oo as can be clearly 
seen from expressions (2.18). (The Fourier transform in this case can be regarded in the sense 
of Plancherel.) Therefore, u 11 G L 2 (R 3 ). Clearly, both u 1 and u 11 depend on the choice of the 
partition of unity 1 = go + he., on the choice of R 0 . Since u 11 = Q~ l gof € L 2 {R 3 ) for any J? 0 , 


then 


Q 9of 


u 


II 


L 2 (R 3 ) 


0 as R 0 


+oo. 


— > 0 as Ro — - Too and consequently, 

Il 2 (R 3 ) 

This concludes the proof of proposition 2.3. 

Proposition 2.4. Let f be the Fourier transformation of f on R 3 and f G Li(R 3 ). Then, 
equation (2.16), (2.1b) has a continuous solution u on R 3 that meets boundary condition (2.2). 
The statement of proposition 2.4 is obvious as in this case Q -1 / G Li(R 3 ). 


In our further constructions, however, we will not always be able to guarantee the inclusion 
f G Li(R 3 ) as required in proposition 2.4. Typically, supp f C Di n and we may also assume that f 
is sufficiently smooth on Di n . On the other hand, we do not generally require the differentiability 
of / on the entire R 3 ; f and its derivatives may have the discontinuities (of the first kind) on 
the surface T = dD { n . For any RHS f of this type, we will make sure that when we successively 
approximate f by the smooth functions f n the corresponding smooth solutions u n in a certain 
sense converge to the solution u guaranteed by proposition 2.3. 

Consider a sequence / n , n = 1,2 ,.. ., of infinitely smooth compactly supported on Dj n functions 
that converges to f in the sense of L 2 (D in ): || f - fn\\ L2 ( Diri ) = \\f ~ fn\\ L2 (R^) — ^ 0 as n — > Too. 
(The sequence f n always exists because f G L 2 (Di n ) and the space V(Di n ) of all compactly 
supported infinitely smooth functions on D{ n is everywhere dense in L 2 (Di n ), see, e.g., [43].) The 
Fourier transformation f n of any f n G V( Di n ) is infinitely smooth on R 3 and decays at infinity faster 
than any power of r _1 with all its derivatives. (Fourier transform in the sense of Plancherel obviously 
coincides here with the standard transform in the sense of L\.) Therefore, for any polynomial 
P = P(£, 77, Q: PQ~ l f n G L\{R 3 ) and consequently, the solution u n to the system Lu n = f n is 
also infinitely smooth on R 3 , u n G C°°(R 3 ), and satisfies boundary condition (2.2). We now consider 
the same partition of unity 1 = g 0 + g^ as used when proving proposition 2.3. u n = (q~ x A)" = 

u i + u n-, where u l = ( Q^gohy and ul 1 = (Q _1 Po/») V ; clearly, both u £, G C°°{R 3 ) 
and both u T n , u^ 1 — , 0 when r — - oo. As f n /, then f n 1 f and g 0 f n L2 L^ o) 

g 0 f. Consequently, g 0 f n g 0 f and therefore g 0 f n L — ’ g 0 f. Since Q~ x G L[ oc (R 3 ), 

then Q~ 1 gof n Q~ 1 gof and also for any polynomial P = P(£, rj, (): PQ^ 1 go f n / -— ' 

PQ~ x gof - Therefore, u ^ uniformly converges to u 1 — ( Q ' 1 gofy on R 3 with all its derivatives, 
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q /3“f* t £^ck~\*I3“ f~ ‘’Y t 

d^a^dTi u n dx ^ dyi^dTi u as n — * +°° ( see > e -S-> [44]). As for the second term, obviously 

9ofn } 96 f and consequently, Q^gah L '~^ ) Q~ l gaf- Therefore, u^ 1 L ~ ' u 11 as n — » +oo. 
Thus, we have justified 

Proposition 2.5. Let supp / C D in and f £ L 2 (D in ). For any f n £ V{D in ), the solution 
u n to the system Lu n = f n (see (2.1b)) satisfies boundary conditions (2.2) and u n £ C°°(R 3 ). 


Moreover, if the sequence f n £ V(Di n ), n = 1,2, . . ., converges to f in the sense of L 2 , f n ^ /, 
then each solution u n can be represented as a sum of two terms, u n = u( -f , where u 1 , u n £ 


gg + /?+ 7 


3Q-+/3+7 


C°°(i? 3 ), ui, ul 1 ► 0 as r ► oo, d f a&y0dzl ^ n dxadyPdt , 
+oo. Here u 1 and u 11 are the same as in proposition 2.3. 


u 1 and ui 1 ^ u 11 


as n 


2.3.2. Finite-Domain AP. To implement the DPM-based ABC’s in practice (Section 3), we 
will need to be able to actually calculate the solution to the auxiliary problem. Since the formulation 
of the AP from Section 2.3.1 still involves infinite domain, we replace it by the approximate finite- 
domain formulation that allows easy numerical solution. 

As any linear system of PDE’s with constant coefficients, (2.16) admits the separation of 
variables in the Cartesian coordinates. Therefore, we implement the Fourier transform in the 
cross-stream and span-wise directions: 


(2.20a) u( X ,V,0 =-^ JJ u{x,y,z)e-^y-^dydz 


(2.20b) f(x, V , C) = ±- JJ f(x,y,z)e- i ^- i < z dydz 

— OO 

and obtain a family of one-dimensional systems 

(2.21) C~ + (i V D + i(E - p 2 F - ( 2 H - t q(j) u = f 

that we consider on the entire fine -oo < x < oo for all ( 1 7, Q £ R 2 . Each system (2.21) is 
supplemented by the boundary condition 


(2.22) | ft (a;, 77, £)| < const for - 00 < x < 00, 

which actually implies |tt(ar, 77, ()l — * 0 as \x\ — (compare to (2.2)) if there are no zeros 
among the eigenvalues of the matrix Q = C -1 (ipD + i(E - 77 2 F - ( 2 H - rj(J). The only special 
case, for which the decay of u(x , 77, () when | x | — - 00 cannot be guaranteed, is (77, () = (0, 0); 
therefore, we generally set the boundary conditions in the form (2.22). It, however, has been shown 
in [13] that after the inverse Fourier transform the solution u will still vanish as |.ij increases. 

Note that although designated by the same notations, it and f in formulas (2.20), as well as 
Q, are not the same here as in the previous section. Indeed, the direction x has been left out when 
calculating the Fourier transformations (2.20). This has been done because the natural spatial 
anisotropy prescribed by the direction of the free stream exists in our model and therefore the 
stream-wise coordinate x will be given a special treatment in the finite-domain AP. 



BOUNDARY CONDITIONS FOR 3D EXTERNAL AERODYNAMICS 


17 


Generally, the solution ii(x, 77, Q to problem (2.21), (2.22) can be found as a convolution 


OO 

(2.23) u(x, 77, 0 = j Gi{x - x rj, Qf{x', 77, Qdx' 

— OO 

with the corresponding one-dimensional fundamental solution G\(x, p, (). Then, the solution u to 
Lu — f can be restored by means of the inverse Fourier transform, which eventually yields: 


00 OO 

(2.24) u(x,y,z) = ^ // f G^x-x', 7?, C) 

— OO “OO 

OO 

JJ f(x, y\ z')e- ir >y l - i ^dy'dz'dx , dr 1 d(. 

— OO 

Now, let us consider the new formulation of the AP that is periodic in both cross-stream and 
span-wise directions . Specifically, we introduce the periods Y and Z for the coordinates y and 
respectively, replace the Fourier integrals by the Fourier series, and instead of (2.24) obtain 


(2.25) u YZ (x, y , z) 


V* E E 



X 


/ 

? 


2 irk y 27r k z \ 


OO 

(' P .2nk v , .'2irk v / 

/ / f(x , t/', z')e~ l ~y^ y dy' dz' dx' . 

— OO 


In our previous work (see [13, 17, 18]), we have analyzed the similar periodic formulations for the 
two-dimensional case. It has been shown that for any fixed-size subdomain the periodic solution 
will converge to the original nonperiodic solution as the period increases. These results can be 
transferred to the case of three space dimensions without changes. Namely, let To and Zq be fixed. 
Then, 


(2.26) 


uyz(x, y, z) 


u(x, y, z) as (Y, Z) — » (+oo, +oo), 
, Y 0 Fo z o ^ ^ Zo 

when -- <y< -, __ <2<y . 


The convergence considered in (2.26) is typically uniform. The convergence for the derivatives can 
also be established under some additional conditions (see [18]). Finally, as we are going to solve the 
AP by a finite-difference method (Section 3), certain relations between the period(s) and the grid 
size(s) should hold, see [13, 17] for more detail. We also note that the convergence on a fixed-size 
domain is sufficient for our purposes because for constructing the ABC’s we will need to know the 
solution of the AP only on some neighborhood of the artificial boundary. 

Thus, we have replaced the original infinite-domain AP by the new problem formulated on the 

beam-shaped domain [-oo, oo] X [— Y/2, Y/2] x [~Z/ 2, Z/ 2]. This domain is still infinite in the 

stream- wise direction. To make the entire formulation truly finite, we first introduce some interval 
[0, A"] so that [0, A"] x [— y/2, Y/2] x [-Z/2, Z/2] D Ti. Consequently, systems (2.21) will be 

homogeneous outside [0, A"] for all (pk y ,Ck z ) = Then, boundary condition 
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(2.27a) H (Q(fe) - A(fc)I) u(0,k) = o, 

_ m(k)<o 

prohibits the lion-decreasing modes in the solution of the corresponding homogeneous system as 
x — » —oo and boundary condition 


(2.27b) JJ (Q(k) — A(fc)J) u(X, k) = o, 

5KA(fc)>0 

prohibits the modes that increase as x — > +oo. Therefore, boundary conditions (2.27) are equiv- 
alent to (2.22) in the sense that the solution to (2.21), (2.27) on [0, X ] will be the same as 
the corresponding fragment of the solution given by (2.23). In formulas (2.27), k = (k y , k z ), 
u(-,k) = u(-,r)k y , (k z ), Q(k) i C~ x (ir]k y D + i( kz E - r)l y F - (l z H - i]k y (k z j), A(fc) are the 
eigenvalues of Q(fc), and I is the identity matrix of the appropriate dimension. 

The formulation of the finite-domain AP is therefore complete. For a given compactly 
supported RHS /, supp f C _D* n , it consists of solving system (2.16) on the parallelepiped 
[0, X] X [— Y/2, Y/2] X [—Zf 2, Z/2] with the periodicity boundary conditions in the y and z di- 
rections and boundary conditions (2.27) in the x direction. As mentioned above, by increasing the 
periods Y and Z one can make the solution to this AP arbitrarily close to the original nonperiodic 
solution on any finite fixed neighborhood of D{ n , 

We will designate the Green’s, i.e., inverse, operator of the finite-domain AP by G so that if 
Lu = f then u — Gf. We also introduce the space T 9 f of the RHS’s for the finite-domain AP 
(V/ : supp f C Din) and the space U 3 u of its solutions so that L : U \ — - T and G : T \ — - U . 
Keeping in mind that the functions u 6 U approximate the solutions to the infinite-domain AP of 
Section 2.3.1 in the sense mentioned above, we will henceforth consider those u E U as satisfying 
the appropriate boundary conditions at infinity. 

2.3.3. Generalized Potentials and Boundary Projections. Let us now introduce the 
space of clear traces E. The elements £ £ E are the vector-functions defined on the artificial 
boundary T; typically, for any u 6 U we may consider £ = (u, | ^ where n is the normal to T. 
The concept of clear trace is delineated in [6, 7]. The operator Tr : U \ — r E that associates the 
clear trace with each u £ U is called the clear trace operator. 

Let now some £ £ if be prescribed. One can always find a compactly supported function v 
such that Trv = £. Then, the truncated function f = (Lv)\ D G T can be a RHS for the 
finite-domain AP. The corresponding solution of the finite-domain AP considered only on D ex is 

called the generalized potential with the density £: P£ d = G . The generalized 

potential can be shown to depend only on its density £ and not on the choice of v (see [6, 7]). 

The composition of operators P and Tr, Pp d = Tr P, maps the space of clear traces onto 
itself, Pp : E i — * E. This new operator is a projection, Pp = Pp, and is called the generalized 
boundary projection. Those and only those £ G E that belong to the image of the generalized 
boundary projection, £ G ImPp, or in other words, satisfy the boundary equation with projection 
£ = Pr£, are actually the traces of some u £ U. 

In the next section, we construct the finite- difference counterparts to the generalized potentials 
and boundary projections and apply those to setting the ABC’s. 
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2.4. Computation of the DPM-based ABC’s. 

2.4.1. Formulation of the Difference AP. Let us introduce a Cartesian grid on the par- 
allelepiped [0, X] x [-Y/ 2, F/2] x [0, Z/2] D IV By virtue of symmetry (see Section 2.2 and, in 
particular, formulas (2.14), (2.15) ), we may consider only a half of the domain along the coordinate 
2. The z-grid is uniform with the size h x : x m = mh x , m = 0,1,..., M, Xo — 0, xm ~ X. The 
grids in y and £ can also be uniform with the sizes h y and h z , respectively: yj y ~ —Y/2 + j y h y , 
j y — 0,1,..., 2 J y + 1, yo — —1 /2, y 2 J y +\ — 5 /2, and Zj z — —Zj2 + j z h z , j z = J z ^ . . . , 2 J z + 1, 
zj 2 = —h z / 2, ^2J^+i — Z/2. (For the £-grid, we use here the same indexing of nodes as if it would 
be if we considered the entire interval [—Zj 2, Z/2] rather than only its half [0, Z/2]. This is done 
mostly to keep consistency in the notations.) However, as we expect to have better accuracy for 
bigger periods Y and Z (see Section 2.3.2), it may be convenient for applications to keep the y- 
and 2 -grids uniform only in the vicinity of Di n and then stretch them away from the computational 
domain. This will allow us to cover bigger periods with the same number of nodes. In so doing, we 
can retain the same indexing for the nodes yi y and Zj z but the grid sizes h y and h z will no longer 
be constant. In all our computations (Section 3), we have actually used the stretched grids in the 
y and £ directions. 

We designate the entire three-dimensional Cartesian grid 

by N°, N° = {(x m ,y jy ,z jz )\m = 0,1,..., M, j y = 0, 1, . . ., 2J y + 1, j z = J z , . . . , 2J Z + 1}. The 
solutions u h £ U h of the difference AP will be defined on this grid. We also introduce an- 
other Cartesian grid M°, on which we will define the RHS’s f h £ T h of the difference AP. 
Compared to the nodes of the grid N°, the nodes of the new grid M° are shifted half-size 
in x. M = ^ (x m _ i/ 2 ? yjy •) Zj z ) — 1, . . . , Af, j y — 0, 1, ... , 2 J y T 1, j z — J z , . . . , 2 J z T 1^, where 

1/2 = (m- 1/2 )h x . 

We discretize the operator L of (2.1) on the grid N° with the second order of accuracy. The 
finite-difference scheme is centered with respect to the nodes (m - 1/2 ,j yy j z ). To discretize 
we use the first-order differences in x , this ensures the second order of approximation because the 
residuals are evaluated on the same semi-integer grid M°, on which the RHS’s are specified. For 
the first derivatives and we use the three-point second-order discretization and designate 
the corresponding grid operators by D y and D,, respectively. The dimension of this operators 
is the same as the dimension of the grid because they act on vector-functions u. } j yi . and 
componentwise. On the uniform grid, this discretization turns into the standard central differencing 
as the central node drops out, but if the grid is stretched the discretization contains all three non- 
zero coefficients. The second derivatives |4f, and are discretized by the appropriate 
compositions of the first difference derivatives; D 2 y , Dl, and D y D z , respectively. We will designate 
the discrete direct operator by L h . 

Let now u h = u m j y j z and f h = f m -i/ 2 j y ,j z ' Because of the periodicity in y, 

(2.28a) ti. )0l . = f , o } * = /• J 2 J y +T- 

Also, because of the symmetry/antisymmetry with respect to z — 0 (see boundary conditions 
(2.14), (2.15)) and periodicity in 2 , 


(2.28b) 


tt vJ, = 


“ f -,',2J z +l-j z 5 ~ 0,1,..., J z ; 

/v,A = “/v,2J z +l-j~ ? jz = 0, 1, . . . , Jz, 


for i / 4, 
for i = 4. 


Again, whereas we formally enumerate the ^-nodes from 0 to 2 J z + 1 in (2.28b), this formulas in 
fact show how one can consider only half of these nodes instead. 
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To replace the continuous transforms (2.20), we introduce the discrete transforms T' y> and 

T (z) 

so that for each i : 


(2.29a) 



jz=2J z -\-l jy-2Jy + l 

E E 


Jz — 1 jy— 1 


rp( z ) rpiv) 

kzjzkyjy 


',3y<,3z 


•> 


(2.29b) 


j z =2J z + l jy=2jy-\-l 

fi _ V T^ Z \ ■ P 

J’,ky,k Z Z-^d k z ,j z kyijyJ^jy^jz' 

jz = 1 iy — 1 


The operators T ^ and have the inverse that we denote by T ^ ~ T ^ 1 and = T W 
respectively, so that 


(2.30a) 


kz—Jz ky—Jy 

A — V A T’W T'bd 7V 

‘Jy ijz Z—J Z—d jz,k z jyyky 'ikyik z ? 

k Z — — J Z ky — — Jy 


(2.30b) 


k z =-J Z ky—Jy 


fi _ V'"' V"' /f^) T 1 ^) p 

jyijz Z Z J jz ,k z jyikyJ’jkyikz' 


We require that the operator diagonalize the first and consequently, the second difference 
derivative with respect to y: 


(2.31a) TtoD y TM = diag {i Vky } , f^D 2 y T^ = diag {-^} , 

where i] ky , /c y = — ,/ y , . . . , J y , are real. Similarly, we require that the operator TW diagonalize the 
first and the second difference derivatives with respect to z: 


(2.31b) tWd x tW = diag RU , fWDjrW = diag{-&} , 

where R, k z = — J 2 , . . ., J z , are also real. From (2.31a) and (2.31b) it follows that 

(2.31c) f^T^D y D z T^T^ = -diag {i lky } diagjc*,}. 

Clearly, the columns of the matrix should therefore be the eigenvectors of D y and analogously, 
the columns of TW should be the eigenvectors of D z . 

Note, in practical computations on the stretched grids (Section 3) the eigenvectors and eigen- 
values of D y and D z are calculated with the standard IMSL subroutines. Although the resulting 
bases are, generally speaking, not orthogonal, the accuracy provided by this calculations is high. 
In fact, this accuracy far exceeds any requirements to the accuracy of ABC’s that may originate 
from the accuracy of the interior solver. The inverse operators and 

are also found with the help of the standard IMSL subroutines. 

If, in particular, the grids in y and z are uniform, then T ^ and are reduced to the 
well-known discrete Fourier transforms (from here on, the overbar “ means complex conjugate): 
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(2.32a) 


T^y) 

ky Jy 


g — ikyjyhy ~ 

\/2 Jy + 1 V^J Z + 1 


rp(z) 

kzjz 


g—ikzjzhz 

\f2Jy + \\f2J Z + 1 


(2.32b) 


rp(y) _ _ 

jy,ky ky,jy 


gikyjyhy ^p?* 

\j2Jy + \ \J2J Z + 1 


rp( Z ) 'Jp( Z ) 

jz Mz kzj z 


gik z jz h z 

\j2Jy + 1\/2J Z + 1 


Let us now consider a special class of grids, namely, yj y — ~y 2 j y +\~ 3y for j y = 0,..., J y and 
Zj z = -z 2 j z +i-j z for j z = 0,...,J^, Obviously, all uniform grids belong to this class, for the 
stretched grids it means symmetric stretching. Then, one can make sure that 


(2.33a) 

rp(y) _ 

— ky yjy 

ji(y) 

ky,jy 

li 

fH 

£ 

0 ,.. 



rp( z ) _ 

“ k z ,j z 

f p( z ) 

1 k z ,jz 

II 

M 

f— l 
£ 

0 ,.. 

*5 Jzj 

(2.33b) 

rp(y) 

ky ,2jy-\-l~jy 

= f £, f °r * 

= 0, 

• * * ? Jy ? 


rp(z) 

k z ,2 J z -\-\ — jz 

, 

ii 

s for Jz 

= 0, 

' ' ">J Zl 

and also 






(2.34a) 

jp(y) _ 

jy ; ky 

j'(y) 

3y iky 

II 

0 ,.. 



(2.34b) 


spiv) _ rp{y) 

2Jy + l—jy,ky jy ,ky 

rp( z ) _ rp( Z ) 

1 2J z +l-j Zj k z 1 j z ,k 3 


fol jy 0 , . . . , J y , 
for j z = 0, 


For the discrete Fourier transform on uniform grids, relations (2.33) and (2.34) immediately fol- 
low from (2.32a) and (2.32b), respectively; for the nonuniform grids these relations are verified 
experimentally. 

Substituting (2.33) into (2.29b), taking into account relations (2.28) and also that f h is real, 
we obtain for i ^ 4 (5? means the real part): 


< 2 - 35 »>4„i* 


jy— 2 Jy +1 

jz — 2 Jz 

- 

£ 

V f2T, ( . 9) , , sirW . f.. , ) 

/ ^ \ i^y | Jy |"'z|jz* / Jy J~ / 


iy — 1 

_iz= Jz + l 

- 


1 = -LlM.IM’ = f- 



and for i = 4 (9 means the imaginary part): 


jy — 2Jy + l j z —2J Z 

( 2 - 35 b)/V. i ,|, |M = i •£ E 2 r £’u aT l 

3y—^ jz—Jz+l 

f-,\ky\-\k;\ = ~f-, \ky\,\k Z \’ f--\ky\,\k Z \ = ~ F.fo \,\k z | » /' - \ky | , - \k Z \ = f ■ , \ky | , |*. | ‘ 
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The relations similar to (2.35) can also be obtained for u h on the basis of formulas (2.28), (2.29a), 
and (2.33). Furthermore, taking into account that u h is real and using formulas (2.30a) and (2.34) 
we obtain for the inverse transform if i ^ 4: 


(2.36a) 


JyiJz 


4 "’if £ - **&.*t]X**A *.) + 


kz 1 ky — 1 

h Z —J z 


2 £ 


k Z = 1 


rx,y — <Jy 

2 Y. (^KImU.0 - + 


ky — 1 


171(2) 'T’(y) f.i 

1 jz, O 1 jy, 0 U ; 0,0 


4: 


<i„i. = -* xf x? + 0® 

/f 2 — 1 Ajj; — 1 

l m k k ) 

x z 3 y A y ^y,h z J 


2 if 
&*=1 


Usage of the transforms (2.35) and (2.36) instead of (2.29) and (2.30), respectively, allows us to 
calculate only one fourth out of the total number of coefficients, namely, those for k y = 0, 1, ... , J y 
and k z = 0,1,..., Jz* This obviously implies a four-fold speedup and four-fold shrinkage of the 
storage requirements when implementing in practice the separation of variables for the difference 
AP. 

In the transformed space, instead of L h u h = f h we obtain a family of one- dimensional systems: 


(2.37) 


where 


Akilm,k + •Bfc’W'm— l,fc “ fm-l/2,ki 

m — 1, . . . , M, k = (& y , fc 2 ), 
ky — 0, . . . , Jy, k z = 0, . . . , J z , 


(2.38) 


Ak -K C+ 2 


tr !ky „ , *0fe 


£> + 


171 171 I T ^7ky Ck 2 , 

2 


B k = -—C + 
h T 


ir lky r, , *C fc 


-£> + 


Vk,, „ Cl TT VkyCkz 


2 2 2 


For each system (2.37), we have to specify the boundary conditions at m = 0 and m = M. 
Analogously, to the continuous boundary conditions (2.27), the boundary conditions for the discrete 
system should explicitly prohibit the corresponding growing modes of the solution. This can be 
achieved by setting 
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(2.39a) 


and 


n 

|A(fc)|>l 


(Qk - A(ib)J) 


uo,k = O, 


(2.39b) 


n 

|A(fc)|<l 


(Qk - X(k)I) 


UM,k = O, 


where Qk = A k 1 Bk , \(k) are the eigenvalues of Qk - and I is the identity matrix of the appropriate 
dimension. 

The finite- difference AP has thus been formulated completely. It consists of solving the discrete 
system L h u h = f h on the grid N° with the RHS specified on the grid M°. The boundary conditions 
in the directions y and 2 are periodicity and symmetry, see (2.28). The boundary conditions in the 
direction x are specified by formulas (2.39) in the transformed space separately for each component 
after the original system L h u h = f h has been reduced to (2.37), (2.38) by the separation of variables 
(2.35), (2.36). The methodology for solving systems (2.37), (2.38) with boundary conditions (2.39), 
as well as the specific structure of these boundary conditions, are studied in the next section. 

2.4.2. Solvability of the Difference AP. Let us first concentrate here on the incompressible 
case, when the 4x4 system matrices are given in (2.1b). For simplicity, we will temporarily omit 
the indices k. If r] = rjk y / 0 or ( e / 0, then the solutions A 5 = A 5 (fc) and e s = e 5 (fc), 
s = 1, . . .,4, of the problem Bj~e — A A& = o are given by 


(2.40) 


Ar = 


e 4 


1 ?? 2 + C 2 \ / 1 ?? 2 + C 2 

h T 2 Re J \ h T 2 Re 


A;< = 


A 4 = 


\A? 2 + C 2 _ ( vV + C 2 , _i 

2 h x ) \ 2 K 


\hFTc + i ) (yftF+e 1 


-l 

= a 2 

r 

. -1 


2 


Ilr 


2 


hr 


Cl = [o, o, -C, vf 


e 2 


«3 = 


0 , 1 , 


— ti) 


1 1 


Re ’ i?e 


-y/v 2 + C 2 + - Re " ) -iv> ~*'C 

y v 2 + ( 2 + r/ , -yV + C 2 , -iv, -*( 


1 1 


From (2.40) we see that we have to analyze two different cases. In the regular case when 
\A? 2 + C 2 /2 — 1 /h x 0, none of the eigenvalues A 5 degenerate, the inverse A^ 1 exists, and the 
eigenvalues/eigenvectors (2.40) are also the eigenvalues/eigenvectors of Q The determinant of 
the Gram matrix constructed on the normalized eigenvectors e s from (2.40) can be shown to be 
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(2.41) 


Deto = 




izi±4f 

fie 2 


» 2 +c 2 

fie 2 


1 + 


»? 2 +c 2 

fie 2 



Therefore, if (i] 2 + ( 2 ) / f?e 2 then the eigenvectors e s are linearly independent and for \f rf + £ 2 /2 — 
1/ha; / Owe can diagonalize the system (2.37): 


(2.42) 


Qk&k — diag {A s } , 


where 



^2 63 e i 

e 2 |’ |e 3 |’ |e 4 |. ’ 


Let us note that since rf 2 ^ 1/h 2 and ( 2 ^ 1/h 2 (h y and h z are the smallest grid sizes) then the 

condition (r/ 2 -|- C 2 ) / Re 2 appears not too restrictive. For example, the asymptotic width of 
the plane viscous wake in the far field behind the body is ~ 1/y/Re (see, e.g., [45]). Therefore, 
to resolve this structure it is sufficient to have the grid sizes of order 1 /\[Re as well (Re is an 
effective turbulent Reynolds number), which puts the operator S & of (2.42) far away of the possible 
singularity. We also note that in the formal inviscid limit 1/Re — » 0, the determinant Detcj of the 
Gram matrix (see (2.41)) becomes fully independent of the wavenumbers p and £, which essentially 
means that the “extent of skewness” for the basis {e s } will be constant. 


The solution to the diagonalized system (2.37) is easy to find by marching those components, 
for which X s < 1. from left to right and those components, for which X s > 1, from right to left. It 
is also easy to make sure that boundary conditions (2.39a) essentially imply that the components, 
for which X s > 1, are not specified (i.e., can be arbitrary) at the left end of the interval and the 
components, for which X s < 1 (those that would not decay as m — - — oo), are zero at m = 0. 
Similarly, boundary conditions (2.39b) mean that the components, for which X s < 1, are not 
specified (i.e., any value is admitted) at the right end of the interval and the components, for which 
A s > 1 (those that would increase as m — * -foo), are zero at m = M. Let us also note that 
1/Re may be arbitrarily small but as long as it is positive, |A S | ^ 1 for all s. Consequently, we 
have only growing and decaying modes and no constant or oscillating modes in the solution of 
the corresponding homogeneous system. Therefore, in accordance with the results of [46] we have 
arrived at 


Proposition 2.6. Let r/k y / 0 or Ck z i 1 0; let also + Cfc 2 / 2 - l/h x ^ 0. Then, 
system (2.37), (2.38), (2.1b) with boundary conditions (2.39) is uniquely solvable and well- 
posed for any compactly supported RHS f m -i/ 2 ,k- The constant in the well-posedness estimate 


In, fc || < const 


f k does not depend on M . 


Note, the system (2.37), (2.38), (2.1b), (2.39) can also be solved using the methodology of [16]. 


The case \/tj 2 + (" 2 /2 — 1 /h x = 0 requires special analysis. In this case, A3 = 0 and also formally 
A4 = 00. In fact, however, it is easy to make sure that both matrices A), and are singular for 
\A? 2 + C 2 /2 — l/hx — 0. Let us therefore consider a regular pencil of matrices Ak + pBk (see, e.g., 
[47]). We can rewrite these pencil as follows: Ak+pB^ - (Ak - Bk) + (p+l)Bj. = A' k + (p + l)Bk- 
As A' k — , this matrix is non-singular and therefore Ak + pBk = A' k + (p + 1 )A' k ~ 1 Bkj- 

The combination of matrices in the brackets can be diagonalized, which yields: 
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(2.43) A k + nB k = 


A' k S' k 


1 

0 

0 

0 


0 

(? + h^Re) + A* (-5 + hTm) 
0 
0 


(3 + h ^ m ) + /* (~3 + 

0 


1 ' 

h x Re j 


0 
0 
0 

“A* J 


Si, 


-1 


where S' k is the corresponding similarity transform. (It is easy to make sure that all the eigenvectors 
are linearly independent so that non-singular S k does exist.) 

From representation (2.43) we conclude that there are still three components in the solution that 
should be calculated by marching from left to right and one component that should be calculated 
by marching from right to left. This obviously matches the structure of boundary conditions (2.39) 
as the latter can also be multiplied from the left by a nonsingular matrix A k . In fact, the pencil 
A k + )iB k has one zero elementary divisor that corresponds to marching from right to left, at least 
one “infinite” elementary divisor that corresponds to marching from left to right, and may have 
either two finite elementary divisors or another two “infinite” elementary divisors that would also 
correspond to marching from left to right. Clearly, any of these marching procedures will easily 
lead to an M-independent estimate of the resulting solution via the prescribed RHS. Therefore, we 
have justified 


Proposition 2.7. Let r) ky ^ 0 or ( kz ± 0; let also yjrg y + Q,J2 - l/h x = 0. Then, 
system (2.37), (2.38), (2.1b) with boundary conditions (2.39) is uniquely solvable and well- 
posed for any compactly supported RHS f m -i/2 t k> The constant in the well-posedness estimate 


||tt. ? fc|| < const 


f M 


does not depend on M . 


Let us now mention that for the discrete Fourier transforms on the uniform grids rj ky ~ 
sin j fay and ( kz = sin ^ 27r/ ^^ Then, to avoid the considerations that result in propo- 

sition 2.7 and to restrict oneself by the case of proposition 2.6 only, one can impose the following 
limitation on the grid sizes: hf 2 > (h~ 2 +hf 2 )/ 4. We also note that the general analysis of constant- 
coefficient ordinary difference equations based on the canonical forms of the corresponding pencils 
of matrices can be found in work [48]. 


The analysis of the last remaining case, rj ky = ( ky = 0 k = o, is straightforward as 
Qo — —I and the solution of (2.37), (2.39) can therefore be found by marching all the components 
from left to right. In accordance with the results of [46] and [16], the well-posedness constant in 
this case is proportional to M . 


For the compressible case (2.1c), the similar results also hold. However, the analytical expres- 
sions of type (2.40) are generally hard to obtain, so the actual eigenvalues and eigenvectors must be 
calculated numerically (we again use the standard IMSL subroutines). The critical value, for which 
the eigen-basis becomes singular (see proposition 2.7) is now ^fv ky + ( 2 z /2 — yj 1 - M^/h x = 0. 

We will designate the Green’s, i.e., inverse, operator of the difference AP by G h so that if 
f h € T h and L h u h = f h then u h = G h f h and u h £ U h . 


2*4.3. Difference Potentials and Projections. Let St m-i/2j y ,j z be the stencil of the dif- 
ference operator L h ; according to Section 2.4.1, we use the first-order differences for the coordinate 
x and the central-type differences and their products for the coordinates y and z. Let us also 
introduce the following grid sets (the overbar Di n here means the set-theoretical closure): 
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Fig. 2.2. Continuous artificial boundary T , grid boundary 7, collocation grid a on and ghost nodes Ti for a 
typical three-dimensional configuration. 


(2.44) 
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By definition (2.44), M, n and M ea; do not have common nodes. The sets N m and N er already have 
some common nodes because these sets are swept by the stencil St m _ 1 / 2//!( , ? . as it is applied to every 
node from M rn and M ea ,, respectively. The intersection of N m and N ex is called the grid boundary 
7. It is actually a multi-layered fringe of nodes of the auxiliary Cartesian grid concentrated in the 
vicinity of the continuous artificial boundary T. Similarly to the continuous case (Section 2.3.3), the 
density of the generalized difference potential will be defined on the grid boundary 7. An example 
of the grid boundary (actually, a few planar cross-sections of this set) for a typical configuration 
studied in this paper is shown in Figure 2.2. 

The difference clear traces £ Sy of the functions u h £ U h are now defined as merely the 


= $y, Tr h : U h 


contractions to the grid boundary, i.e., Tr h u h d = u 1 

Let now some £ 7 £ be prescribed and v h be a grid function defined on N° such that 
Tr h v h = £ 7 . Clearly, there are many functions v h that would meet this condition, for example 


v h = f £7 on 7, 

( o on N°\7- 


(2.45) 
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Then, consider the function 


(2.46) J* 3 = { Lhvh 

I O on Me*, 

where v h is defined by (2.45), and solve the difference AP with this RHS f h of (2.46). The resulting 
solution considered only on N ex is called the generalized difference potential with the density £ 7 : 


(2.47) 


d = f 



N 


ex 


(f h in (2.47) is defined by (2.46)). Analogously to the continuous case (Section 2.3.3), the general- 
ized potential P h £ 7 of (2.47) can be shown to depend only on its density £ 7 and not on the choice 
of v h (if the latter differs from (2.45) but still Tr h v h = £ 7 ), see [7]. 

The composition of operators P h and Tr h , P 7 d = Tr k P h , maps the space of the difference 

clear traces onto itself, P^ : Sy ■ — - E 1 . This new operator is a projection, Pi’ 2 = Pp , and is called 
the generalized difference boundary projection. Those and only those £ 7 G Ey that belong to the 
image of the generalized difference boundary projection, £ 7 G ImP 7 , or in other words, satisfy the 
difference boundary equation with projection £ 7 = Pi ’ £ 7 , are actually the traces of some u h G U h . 


Note, numerical verification of the projection property P 7 2 = P^ is an ideal test for accuracy 
of the solution of difference AP. In our practical computations for different geometries on different 


grids, we have always been able to obtain for arbitrary £ 7 ’s: 


- Plli 


7 S7 


< 1(T 9 . This, in 


particular, justifies the usage of the stretched grids when solving the difference AP. 

As mentioned before (Section 2.3.2), we consider the continuous functions u £ U as satisfying 
the appropriate boundary conditions at infinity because the difference between the non-periodic 
solution and its periodic approximation is controlled by Y and Z and can, in fact, be made as 
small as initially prescribed. The discrete space U h approximates the continuous space Z7, therefore 
we consider those grid densities that belong to ImP 7 , £ 7 £ Im Pjf as admitting the exterior 
complement in the right sense. In other words, these and only these functions admit such a 

complement u^ x P h & y that satisfies the boundary conditions of the difference AP (see (2.28), 

(2.39)); this complement can therefore be made arbitrarily close (near D { n ) to an original linearized 
exterior solution; in the next section, it is used for setting the difference ABC’s. 


2.4.4. Global DPM-based Artificial Boundary Conditions. Having constructed the 
procedure for calculating the generalized difference potentials and projections, we can now pro- 
vide for a closure to the discretized Navier-Stokes system that is solved inside the computational 
domain Di n , i.e., obtain the ABC’s. As mentioned in Section 2.2, the interior solvers typically 
involve some sort of pseudo-time iterations. To make every step of the iteration procedure, we 
need to know the previous-step solution everywhere on the grid, including the ghost nodes Ti. If 
these data are available, then on the next step we will know the solution everywhere except on 
T i. Consequently, to advance another iteration we will have to supplement the missing data on T\. 
This will be done by projecting the available boundary data at T onto the ‘Tight manifold”, i.e., 
the one that admits the right exterior complement (see the previous section), and then calculating 
this complement on P 1( In so doing, we can obtain the missing relations between the values of the 
solution on T and Ti every time the ghost nodes need to be updated. 

First, let us introduce the intermediate collocation grids a and oq on both surfaces T and Ti. 
An example of such a C T is shown in Figure 2.2. These grids are typically a few times coarser than 
T and r i. Usage of the collocation grid on T is an element of general procedure of the difference 
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potentials method [7]. Moreover, for the specific problem under study usage of the collocation grids 
results in multi-fold acceleration of the computational procedure and also in memory savings if the 
boundary conditions are implemented in the matrix form (see below). 

Then, let us take £p = (u, where n is the normal to T (these data are available from 

inside the computational domain D{ n on every iteration) and, using the clustering R a on T, obtain 
£a* The latter procedure (clustering), in fact, implies local averaging or smoothening along T. 
Furthermore, we drop normals from all nodes 7 to the surface T and interpolate with sufficiently 
high order to the feet of these normals. The corresponding operation is denoted R n ; typically, we 
use the bi-quadratic surface interpolation. Having obtained u and at the feet of the normals, 
we use the first two terms of the Taylor expansion (denoted 7 r 7 ) and obtain £ 


7* 


(2.48) 


£7 — ^"7 R n R (t ( 


Then, we calculate the potential P h C i for the density = P^ 7 and interpolate it (operation R ai ) 
from N ea: to the nodes ai C fp 


(2.49) 


u 


_ T>h,£t — D E>/l 1 


<71 


^<y 1 


= R ai p n i 


■7’ 


The second equality in (2.49) holds because of the projection property of . Finally, the missing 
values of the solution at the nodes Ti are obtained from u by means of interpolation along the 

G\ 

surface Ti, which altogether yields the nonlocal DPM-based ABC’s in the form 


(2.50) 







here the operation T is composed of the operations (2.48), (2.49), and interpolation along Tj. As 
mentioned above, in the course of the iteration procedure boundary condition (2.50) is applied 
every time we need to update the values of the solution at the ghost nodes T 1. The implementation 
of ABC’s (2.50) can either be direct or involve preliminary calculation of the matrix T. In the 
latter case, the runtime implementation of the ABC’s (2.50) is reduced to a matrix-vector multi- 
plication. Moreover, in this case we can do the first clustering R a and the last interpolation along 
Ti separately, i.e., leave these operations out of the structure of T. Then, instead of (2.50) one can 
write 


u 


(71 


= T% 


where both the dimension of T* and its computational cost are many times smaller than those of 
T from (2.50). 

Let us also note that we need to know the potential only on some neighborhood of the surface 
Fi (see (2.49)). At the same time, according to (2.45) and (2.46) the density of the potential differs 
from zero only near 7. Therefore, for both direct T^ y \ (see (2.29)) and inverse T^ y \ 

(see (2.30)) transforms we actually have to take into account only a few “non-zero” nodes out of 
the total numbers of 2 J y + 1 and J z + 1 along y and z y respectively. This effectively makes the 
computational cost of these transforms to grow linearly rather than quadra.tically with respect to 
2 J y + 1 and J z + 1, and obviously implies a very substantial reduction of the required computer 
resources. 
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3. Numerical Results. 

3.1. Two-Dimensional Summary. For the reason of completeness, we first briefly comment 
on the two-dimensional results from our previous work (see [14, 15, 18]). In that work, we have 
calculated the subsonic and transonic viscous flows past single-element airfoils (NACA0012 and 
RAE2822). 

The two-dimensional computational domain is formed by the C-type curvilinear grid generated 
around the airfoil. On this grid, the Navier- Stokes equations are integrated using the code FLOMG 
by Swanson and Turkel [10, 11, 12]. The standard treatment of the external boundary in the code 
FLOMG is based on the locally one-dimensional characteristics analysis, which may or may not be 
supplemented by the point- vortex correction [32]. 

Basic conclusions that could be drawn from our two-dimensional numerical experience are the 
following. The DPM-based ABC’s are geometrically universal, algorithmically simple and easy to 
implement along with the existing solver. For the large computational domains (30-50 chords of the 
airfoil), the performance of the standard methods and the DPM-based ABC’s is very close to one 
another. As, however, the artificial boundary approaches the airfoil the discrepancy between the 
corresponding solutions increases. The lift and drag coefficients obtained on the basis of the two- 
dimensional version of boundary conditions (2.50) deviate from their asymptotic (50 chords) values 
much slighter (within fractions of one percent) than the coefficients obtained with the local ABC’s 
do. In other words, the nonlocal DPM-based ABC’s allow one to use much smaller computational 
domains (as small as 2-3 chords) than the standard boundary conditions do and to still maintain 
high accuracy of of the numerical solution. Moreover, if we compare three models: DPM-based, 
point- vortex, and standard local (characteristics-based), then it turns out that the DPM-based 
ABC’s display the best performance for small computational domains, the performance of the local 
characteristic boundary conditions for small domains is very poor, and the point-vortex boundary 
conditions perform much better for the lift than they do for the drag coefficient. This behavior 
seems reasonable since the point-vortex model is a lift-based treatment. 

We also note that for certain variants of computation the DPM-based ABC’s may noticeably 
speed up (by up to a factor of three) the convergence of the multigrid iterations, see [13, 14, 15]. 
Some discussion on combined implementation of the DPM-based ABC’s with multigrid is contained 
in [18]. 


3.2. Three-Dimensional Computations. The DPM-based boundary conditions (2.50) 
have been combined with the interior Navier- Stokes solver and used for calculating viscous flows 
around the ONERA M6 wing for different regimes that range from very low to transonic Mach 
numbers and include both attached and separated turbulent flows. 

We use the NASA-developed code TLNS3D by Vatsa, et al. [21] to integrate the thin-layer 
equations on the curvilinear C-0 grid (see Figures 2.1 and 2.2) generated around the wing. The 
code is based on the central-difference finite-volume discretization in space with the first- and 
third-order artificial dissipation. The steady-state solution is obtained by means of a pseudo-time 
iteration procedure; the integration in time is done by the five-stage Runge-Kutta algorithm (with 
the Courant number calculated locally) supplemented by the residual smoothing. For the purpose 
of accelerating the convergence, the multigrid methodology is implemented; in our computations 
we used three subsequent grid levels with V cycles; the full multigrid methodology (FMG) could 
be employed as well. In addition, we use the preconditioning technique of [49] to improve the 
convergence to steady state. We implement the DPM-based ABC’s (2.50) only on the finest level 
of multigrid on the final FMG stage; the boundary data for coarser levels are provided by the 
coarsening procedure. Moreover, even on the finest level we implement the DPM-based ABC’s only 
on the first and the last Runge-Kutta stages, which has been shown to make very little difference 
compared to the implementation on all five stages; the boundary data for the three intermediate 
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stages are provided from the DPM-based ABC’s on the first stage. Unlike the two-dimensional case, 
the standard treatment of the external boundary in three dimensions (code TLNS3D) is based merely 
on the locally one- dimensional characteristics analysis and extrapolation (as the point-vortex model 
is not applicable). 

All three-dimensional flows that we have analyzed are turbulent. In the near field (i.e., inside 
the Na.vier- Stokes solver is supplemented by a special turbulence model to account for the 
corresponding phenomena. Depending on the specific flow variant, either an algebraic or a differen- 
tial turbulence model can be employed. In the far field, we use Boussinesq’s concept of the effective 
turbulent viscosity, i.e., effective Reynolds number (see [15]). This simplest approach has been 
found to produce accurate results when incorporated in the structure of the DPM-based ABC’s. 
The value of the Prandtl number for all the calculations was either Pr = 0.72 (air) or Pr = 1. 

In all the cases below, the auxiliary Cartesian grids are stretched along the coordinates y and z. 
The stretching typically starts outside T i; the stretching factors (we use the geometric progression) 
vary between 1.07 and 1.1 for different variants. The typical values of Y and Z that we have used 
also vary between 20 and 30 sizes of Di n in the cross-stream direction and 4 and 10 sizes of D{ n in 
the span-wise direction. The uniform Cartesian grid in the vicinity of Di n is always chosen so that 
the distance between T and is well resolved. 

We should also emphasize that in spite of their nonlocal nature the DPM-based ABC’s (2.50) 
are geometrically universal. In other words, these boundary conditions can be obtained for the 
boundary T of any irregular shape by means of the same computational procedure. This conclusion 
directly follows from the previous considerations and has also been repeatedly corroborated in the 
numerical experiments. Moreover, ABC’s (2.50) appear easy to incorporate in the structure of the 
existing flow solvers, which has been corroborated in practice as well, and which is very important 
from the standpoint of applications. The issue of the computational cost of boundary conditions 
(2.50) and some possible ways of its reduction will be addressed later on, in Section 3.2.4. 

3.2.1. Low Mach Number Regime. We first consider a very low speed flow, M 0 = 0.01, 
which, in fact, is close to the truly incompressible case. Preconditioning [49] makes the analysis of 
this flow possible with TLNS3D. The flow is turbulent with the molecular Reynolds number based on 
the root chord of the wing Re o = 11.7 * 10 6 ; the angle of attack is a = 3.06°; there is no separation 
and the turbulence inside D{ n is simulated using the Baldwin-Lomax algebraic model, which is 
based on the concept of mixing length. 

Since the free-stream Mach number is so small, we have implemented here the incompressible 
version of the nonlocal DPM-based ABC’s (2.50) constructed on the basis of matrices (2.1b). In 
Table 3.1, we present the results of calculations for two different computational domains of the 
“average radii” of 10 and 1.25 root chords of the wing, respectively (root chord means the chord 
length at = 0). 


Table 3.1 

ONERA M6 ; Af 0 = 0.01; Re 0 = 11.7 • 10 6 ; a = 3.06°. 


“Average radius” of Di n 

1.25 root chords 

10 root chords 

Dimension of the grid 

197 x 49 x 33 

Type of ABC’s 

standard 

DPM 

standard 

DPM 

Full lift, C L 

0.2052 

0.1954 

0.1940 

0.1939 

Relative error 

5.78% 

0.77% 

0% 

0% 

Full drag, Cp X 100 

0.695 

0.685 

0.681 

0.681 

Relative error 

2.1% 

0.58% 

0% 

0% 


In both cases, we used the C-0 type grids of the same dimension 197 X 49 X 33; for the small 
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domain the grid was obtained by scaling down the bigger grid and was obviously finer in the near 
field. One can see that for the big domain the results (force coefficients Cl and Co) obtained using 
both methods are very close to each other. However, as the domain shrinks the accuracy obtained 
with the DPM-based procedure appears much better than the accuracy provided by the standard 
methodology. In other words, the nonlocal DPM-based ABC’s (2.50) allow one to essentially reduce 
the size of the computational domain without compromising the accuracy. This confirms that if 
the structure of the far-field solution is correctly taken into account by means of the ABC’s then 
within a certain range of domain sizes the computed near-field solution becomes essentially domain- 
independent. We also note that as the near- field grid for the small domain is finer than for the big 
domain then the associated truncation error can be expected smaller. 

3.2.2. Subsonic Regime. The next case is a subcritical (i.e., fully subsonic) compressible 
flow for Mo = 0.5. Here, the free-stream Mach number is already high enough to make the 
compressibility effects very essential but on the other hand, it is still not too high so that the flow 
remains locally subsonic throughout the entire domain. The angle of attack and the molecular 
Reynolds number for this case are the same as for the previous one: a = 3.06°, Re 0 = 11.7 • 10 6 . 
The flow is also fully attached and the turbulence model inside D{ n is algebraic (Baldwin-Lomax). 

The DPM-based ABC’s (2.50) for this case were constructed on the basis of non- symmetrized 
matrices (2.1c). For this specific value of Mach number, Mo — 0.5, the “extent of non-symmetry” in 
the system matrices (2.1c) still appears quite acceptable. However, for low Mach numbers Mo < 0.1 

treated in the compressible framework (unlike in Section 3.2.1), the usage of symmetrizer (2.3) and 
matrices (2.4) can be recommended. On the other hand, we should note that in work [18] we 
have been able to obtain accurate two-dimensional results for M 0 = 0.01 without symmetrizing the 
system matrices in boundary conditions. 

In Table 3.2, we compare the results of calculations for three different computational domains. 

Table 3.2 

ONERA M6: Mo = 0.5; Re 0 = 11.7 • 10 6 ; a = 3.06°. 


“Average radius” of 

1.25 root chords 

2 root chords 

10 root chords 

Dimension of the grid 

197 x 49 x 33 

Type of ABC’s 

standard 

DPM 

standard 

DPM 

standard 

DPM 

Full lift, C L 

0.2218 

0.2065 

0.2185 

0.2065 

0.2081 

0.2072 

Relative error 

6.58% 

0.34% 

5.0% 

0.34% 

0% 

0% 

Full drag, Cjj x 100 

0.817 

0.791 

0.793 

0.791 

0.787 

0.788 

Relative error 

3.8% 

0.38% 

0.76% 

0.38% 

0% 

0% 


Like in the previous case (Section 3.2.1), here the DPM-based ABC’s produce much more accu- 
rate solutions on the small computational domains than standard boundary conditions do. This 
obviously amounts to either saving the computer resources while preserving the accuracy of compu- 
tations or improving the accuracy while keeping the computational cost at the same level. Of course, 
lower levels of the truncation error for finer grids on the small domains can also be anticipated here 
as in the foregoing low Mach number case. 

3.2.3. Transonic Regime. Most of the standard test cases for flows around the ONERA M6 
wing are transonic (see, e.g., the experimental work [50]). In such flows the free-stream Mach 
number is sufficiently high so that the local speed exceeds the speed of sound in some bounded 
region near , the upper surface of the wing. This leads to the formation of a supersonic (i.e., 
supercritical) “bubble”, which typically has a sonic-surface type upstream boundary and a shock- 
wave type downstream boundary. 
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Attached Flow. The first transonic case that we present is Mq = 0.84, a — 3.06°, Re o = 
11.7 • 10 6 . In this case, the angle of attack a remains sufficiently small so that the weak shock 
;t the upper surface of the wing does not cause the flow separation. Therefore, we still use the 
Baldwin-Lomax model for simulating the turbulence inside D{ n . An important difference compared 
to the previously studied cases is that here we cannot bring the artificial boundary as close to the 
wing as done in Sections 3.2.1 and 3.2.2. The reason is that our far-field treatment is purely subsonic 
and therefore, the artificial boundary should not come to the boundary of the supercritical bubble 
too close. Therefore, we ran our computations for two domains, the “radius” of the big one is still 
about 10 root chords of the wing and the “radius” of the small one is about 3 root chords of the 
wing. Unlike in the previous cases, here we constructed the C-0 grids of different dimension for the 
domains of different size; the smaller (3 root chords) grid is an exact subset of the bigger (10 root 
chords) grid. This has been done in order to completely eliminate any influence that the change of 
the grid in the near field may possibly exert on the calculated solution. 

The nonlocal ABC’s (2.50) for this case were again constructed on the basis of matrices (2.1c). 
In Table 3.3, we compare the computed results (calculated lift Cl and drag Co coefficients) for 
two different types of ABC’s on two different domains. 


Table 3.3 

ONERA M6: M 0 = 0.84; Re 0 = 11.7 * 10 6 ; a = 3.06°. 


“Average radius” of D\ n 

3 root chords 

10 root chords 

Dimension of the grid 

197 x 49 x 33 

209 x 57 x 33 

Type of ABC’s 

standard 

DPM 

standard 

DPM 

Full lift, C L 

0.298T0.004 

0.2798 

0.2805 

0.2786 

Relative error 

6.24%±1.43% 

0.43% 

0% 

0% 

Full drag, C'd x 10 

0.168T0.008 

0.1537 

0.1542 

0.1531 

Relative error 

8.95%±5.19% 

0.39% 

0% 

0% 


For the small computational domain, the DPM-based ABC’s again clearly outperform the standard 
method from the standpoint of accuracy. We also note that in this case the total number of nodes 
in the bigger grid is about 25% more than in the smaller grid, which obviously implies an accordant 
increase of the associated cost of computations. 

Even more important, for the transonic case the DPM-based ABC’s influence not only the final 
accuracy of the solution but also the convergence rate of the iteration procedure employed inside 
Di n . Namely, for the standard ABC’s the multigrid iterations on the small computational domain 
converge noticeably slower than they do for the DPM-based ABC’s. In fact, for the same 500 V- 
cycles on the finest multigrid level, we simply have not been able to obtain a stable solution for the 
3 root chords domain with the standard boundary conditions. That’s why the corresponding data 
in Table 3.3 are given with the error bands indicated. The convergence history for the transonic 
computations on the 3 chords domain is given in Figure 3.1a. for the residual of the continuity 
equation and in Figure 3.1b for the number of supersonic points in the domain. Note, the latter 
quantity is deemed very sensitive for calculation of the transonic flows. 
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Fig. 3.1a. ONERA M6 : Mo — 0.84, Re o = 11.7 • 10 6 , a = 3.06°. Convergence history for the residual of the 
continuity equation. (( Average radius” of D tn is 3 root chords of the wing; the dimension of the grid is 197 x 49 x 33. 



Work 

Fig. 3.1b. ONERA M6: Mo = 0.84, Re o = 11.7 • 10 6 , a = 3.06°. Convergence history for the number of 
supersonic nodes in the domain. 1 Average radius ” of Di n is 3 root chords of the wing; the dimension of the grid is 
197 x 49 x 33. 

From Figures 3.1 one can easily see that the difference in the multigrid convergence rates for the 
different types of ABC’s can be as big as approximately a factor of three. 

The history of convergence of the same two quantities for the big (10 root chords) computational 
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domain is presented in Figures 3.2. 



Fig. 3.2a. ONERA M6: Mq = 0.84, Re o = 11.7 * 10 e , a = 3.06°. Convergence history for the residual of the 
continuity equation. “ Average radius ” of Di n is 10 root chords of the wing; the dimension of the grid is 209 x 57 x 33. 


& 

c 

0 

01 


o 

'c 

o 

£ 


Q 

Q_ 

3 


C/D 


32000 


28000 


24000 


20000 


16000 



Work 


Fig. 3.2b. ONERA M6: Mo = 0.84, Re o = 11.7 • 10 6 , a = 3.06°. Convergence history for the number of 
supersonic nodes in the domain. Average radius” of Di n is 10 root chords of the wing ; the dimension of the grid is 
209 x 57 x 33. 


We see that that in this case the DPM-based ABC’s also provide for some convergence speedup 
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although the difference between the two methodologies appears less dramatic. This seems reason- 
able because one could generally expect that the bigger the computational domain, the smaller is 
the influence that the external boundary conditions exert on the numerical procedure. 

Let us also note that on the small (3 root chords) domain the two algorithms apparently 
converge to quite different solutions (this is most clearly seen in Figure 3.1b), whereas Figure 3.2b 
allows one to assume that on the big (10 root chords) domain the final solutions are close to one 
another. The data from Table 3.3 corroborate these conclusions. This behavior of the solution again 
fits into the aforementioned concept that the overall impact of the ABC’s on the computational 
algorithm decreases as the domain enlarges. 

Separated Flow. When one increases the angle of attack a in the transonic regime, the flow 
pattern changes. The shock on the upper surface of the wing becomes stronger. Since the chord 
length of the wing decreases span-wise as z increases (see Figure 2.1), then the stream-wise size 
of the supersonic bubble decreases as well, and eventually the upstream sonic surface and the 
downstream shock wave meet somewhere in the area, close to the wingtip. For sufficiently strong 
shocks this, in particular, produces flow separation on the upper surface of the wing. We have 
analyzed the separated flow of this type for Mo = 0.84, a = 5.06°, Re o = 11.7 • 10 6 . 

The separation zone on the upper surface of the wing for this case is relatively small, the flow 
fully re-attaches before the trailing edge so that no phenomena associated with the separation are 
present in the wake. However, the simulation of such flows already requires more sophisticated 
turbulence models inside the computational domain; we have used the the two-equation Menter’s 
model [51]. Moreover, it requires much finer grids in the near field than the simulation of the 
attached flows does. 

As in the previous transonic case, the global ABC’s (2.50) are constructed here on the basis of 
matrices (2.1c). The computations are conducted for two different domains of the “average radii” 
of 3 and 10 root chords of the wing, respectively, on the grids of the same dimension 193 X 49 x 33; 
the smaller grid is obtained by scaling down the bigger grid (analogously to Sections 3.2.1 and 
3.2.2). In Figure 3.3, we present the distribution of the pressure coefficient 


^ _ P-Po 

” 1 9 

\po‘Uo 2 


on the upper and lower surfaces of the wing in the cross-section £ = const at the 90% of semi-span. 
The 90% of semi-span station corresponds to the area of developed separation. The three solutions 
that we have computed in this case are for the global DPM-based ABC’s on the 3 and 10 root 
chords domains and standard ABC’s on the 10 chords domain. These solutions are compared in 
Figure 3.3 against the experimental data. 

From Figure 3.3 we conclude that all three numerical solutions very well match one another 
and also match the experimental data to a reasonable degree of accuracy. We also emphasize 
that analogously to the previous cases, the DPM-based global ABC’s (2.50) are quite capable of 
generating an accurate numerical solution on the small domain for this separated flow case. On 
the other hand, unlike in the previous cases, the standard ABC’s are simply unable to produce 
a convergent solution on the 3 root chords computational domain for the a = 5.06° separated 
flow around the ONERA M6 wing. In other words, the multigrid algorithm with the standard 
ABC’s fails to converge on the small computational domain. Let us recall that for an easier case of 
a — 3.06° the convergence of the standard procedure on the small domain was only slowed down 
but not completely destroyed. The history of convergence of the residual of continuity equation for 
the case a = 5.06° on the small domain is presented in Figure 3.4a. 
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Fig. 3.3. ONER A M6: Mo = 0.84, Re o = 11.7 • 10 6 , a = 5.06°. Surface pressure distribution at the 90% of 
semi-span (x/c: x is the coordinate calculated from the leading edge , c is the local chord length ). Dimension of all 
grids is 193 x 49 x 33. 



Fig. 3.4a. ONERA M6: Mo = 0.84, Re o — 11.7 • 10 6 , a = 5.06°. Convergence history for the residual of the 
continuity equation . “Average radius” of Di n is 3 root chords of the wing; the dimension of the grid is 197 x 49 X 33. 
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At the same time, on the big (10 root chords) domain both algorithms for the a = 5.06° case 
converge at the same rate, see Figure 3.4b. 



Fig. 3.4b. ONERA M6 : Mo — 0.84, Re o = 11.7 * 10 6 } a = 5.06°. Convergence history for the residual of the 
continuity equation. Average radius ” of D tn is 10 root chords of the wing; the dimension of the grid is 197 x 49 X 33, 

Figures 3.3 and 3.4 allow us to conclude that the nonlocal DPM-based ABC’s (2.50) not only 
speed up the convergence of the multigrid iterations but are generally capable of increasing the 
robustness of the entire numerical procedure. 

3.2.4. Computational Cost of the DPM-based ABC’s. In all the three-dimensional 
computations described above, the DPM-based ABC’s were implemented directly, without com- 
puting the matrix of operator T from (2.50). By applying the new procedure only on the first 
and the last R.unge-Kutta stages and only on the finest multigrid level, the total number of the 
required calculations of generalized potential has been brought to a minimum. In so doing, the 
average cost of application of the DPM-based ABC’s (2.50) adds about 20-25% of the CPU time 
to the cost of the same procedure with the standard (characteristics-based) boundary conditions. 
This extra expense is not high (taking into account the improvement of accuracy); moreover, it can 
often be compensated for and even noticeably prevailed over by the convergence acceleration and 
the reduction of the domain size. Besides, to explicitly decrease the computational cost associated 
with the DPM-based ABC’s we plan on the future use of the entry- wise interpolation of boundary 
operators (see [22]) and/or the multiresolution based methodologies (see [20, 22]). We expect that 
the latter can also be employed when implementing the DPM-based ABC’s for multi-block grids. 

4. Conclusions. The new global ABC’s for calculating steady-state external viscous flows in 
three space dimensions have been constructed on the basis of the difference potentials method. The 
approach generalizes and extends our previous two-dimensional results. 

The new ABC’s are capable of greatly reducing the size of the computational domain (compared 
to the standard methods) while still maintaining high accuracy of the numerical solution. This size 
reduction amounts to either the possibility of refining the grid in the near field, which potentially 
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leads to increasing the accuracy, or usage of the smaller-dimension grids while keeping the accuracy 
at the same level. Moreover, the DPM-based ABC’s may noticeably speed up the convergence of 
the multigrid iterations and generally improve the robustness of the entire numerical procedure. 
Finally, the new boundary conditions appear geometrically universal and easy to incorporate in 
the structure of the existing flow solvers. The properties of the new ABC’s have been corroborated 
experimentally by computing the subsonic and transonic flows past the ONERA M6 wing using 
the NASA-developed code TLNS3D. 
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